AD- A 186  127  fllB.flii 


[SECURITY  CLASSIFICATION  or  THIS  PAPE  (When  Pdf.  Entered).  * 

r  REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


EPORT  NUMBER 


a  AFIT/qi/NR  87-84T 

«.  TITLE  (rtnd  Submit} 

On  The  Parabolized  Navler-Stokes  Equations 
Without  Sublayer  Assumptions 


h.  GOVT  ACCESSION  NO.I  »•  RJECIRIENUU:  AT ALOO  NUMBER 

1  U477C  Jl  7 

I.  TYPE  op  REPORT  »  PERIOD  COVERED 

THESIS/DJ2£Jt)7f;;?N 

s.  performing  obg.  report  number 


I 


it.  authoru; 


Stephen  Christian  Pluntze 


3.  CONTRACT  OR  GRANT  NUMBER! a) 


.  PERFORMING  ORGANIZATION  NAME  AND  AOOREtS 

AFIT  STUDENT  AT: 


10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  •  WORK  UNIT  NUMBERS 


Massachusetts  Institute  of  Technology 

I.  CONTROLLING  OFFICE  NAME  AND  AOORESS  ■  <*•  REFORT  DATE 

AFIT/NR  1986 

WPAFB  011  45433-6583  is.  number  of  pages 

_ 147 _ 

.  MONITORING  AGENCY  NAME  a  ADDRESS!!!  olllerenl  from  Central  ling  Ollice)  IS.  SECURITY  CLASS,  (o!  Ihl.t  report) 


UNCLASSIFIED 

Tsa.  DECL  ASSIFICATION/  DO  WNGRADING*' 
SCHEDULE 


.  DISTRIBUTION  STATEMENT  (ol  (M*  Report) 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


17.  DISTRIBUTION  STATEMENT  (ol  the  abstract  an! trod  In  Block  30,  II  dlflerenl  Irom  Raporlj) 


DTIO 

,ELECTE| 
NOV  0  4 19871 


IS.  supplementary  notes 

APPROVED  FOR  PUBLIC  RELEASE:  1AW  AFR  190-1 


[  19.  KEY  WORDS  (Continue  on  ravarae  a  Ida  II  neceeeety  and  Identity  by  block  number) 


20.  ABSTRACT  (Conllnu#  on  ravaraa  alda  II  neceeeery  ltd  Identity  by  block  number) 

ATTACHED 


DO  |  JAIWJ  1^73  EOITION  OF  I  NOV  6S  IS  OBSOLETE 


LYNg7  E .  WOLAV  ER  ^  Ji/f  >0 
Dean  for  Research  ana 
Professional  Development 
AFIT/NR 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (HTien  Dale  Enter  ad) 


Enclosed  Is  the  Master  of  Science  Thesis  for  Capt.  Stephen  C.  Pluntze. 
It  Is  147  pages  long  and  was  submitted  to  the  Dept  of  Aeronautical 
and  Astronautlcal  Engineering  at  the  Massachusetts  Institute  of 
Technology  for  the  degree  of  Master  of  Science  in  Aeronautical  and 
Astronautlcal  Engineering. 


Acc«;4on  i  or 


NTIS  CttA&l 
DUG  TAB 


n 


ON  THE  PARABOLIZED  NAVIER  -  STOKES  EQUATIONS 
WITHOUT  SUBLAYER  ASSUMPTIONS 


by 

STEPHEN  CHRISTIAN  PLUNTCE 
Captain,  United  States  Air  Force 

BS,  UNITED  STATES  AIR  FORCE  ACADEMY 
U9£2> 


SUBMITTED  IN  PARTIAL  FULFILLMENT 
OF  THE  REQUIREMENTS  FOR  THE 
*  DEGREE  OF 

MASTER  OF  SCIENCE  IN 
AERONAUTICS  AND  ASTRONAUTICS 


at  tha 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
January,  1987 

©  Stephen  C.  Pluntze,  1987 


The  author  hereby  grants  to  M. I.T.  peraisslon  to  reproduce  and 
to  distribute  copies  of  this  thesis  document  In  vhole  or  in  part. 


Signature  of  Author  _ _ _ 

Department  of  Aeronautics  and  Astronautics 

January  23,  1987 


Certified  by 


Prof.  Judson  R.  Baron,  Thesis  Supervisor 
Department  of  Aeronautics  and  Astronautics! 


Accepted  by _ _ 

Prof.  Harold  Y.  Wachman 
Chairman,  Department  Graduate  Committee 


1 


ON  THE  PARABOLIZED  NAVIER-STOKES  EQUATIONS 
WITHOUT  SUBLAYER  ASSUNPTIONS 

by 

STEPHEN  CHRISTIAN  PLUNTZE 
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ABSTRACT 

''"^A  new  iaplioit.  Iterative  aethod  of  solving  the 
Parabolised  Clavier -Stokes  (PNS)  Equations  olalaa  to  overeoae  the 
elliptic  oharaoter  of  the  eabedded  subaonlo  sublayer  by 
explicitly  introducing  pressure  as  an  additional  state  variable. 

The  Bhutta-Lewis  approach  aakes  no  sublayer  pressure  assumptions. 

The  validity  and  basis  of  that  aethod  is  explored  in  this  thesis 
by  examining  the  relevant  eigenvalues  governing  marching 
stability.  An  original  code  was  also  developed  in  order  to 
examine  the  numerical  character  of  the  marching.  Iterative 
solutions  as  they  develop.  Test  cases  mere  carried  out  for  a  two 
dimensional  wedge  configuration  at  Hach  numbers  3  and  13  and 
Reynolds  numbers  ranging  from  4xlC*  to  ( 1x10*  at.  the  initial  data 
plane.  ifopo  -i0  fc> 

An  eigenvalue  analysis  disclosed  that  the  method  is 
unstable  in  subsonic  regions.’  Introducing  the  additional  state 
variable  does  not  change  the  character  of  the  equations. 

Results  for  the  test  cases  confirmed  the  presence  of 
instability.  Classic  departure  behavior  was  produced  in  tightly 
clustered  grids  and  convergence  to  separated  flow  was  shown  in 
less  clustered  grids.  Marching  was  achieved  only  in  relatively 
high  Reynolds  number  flow  with  a  large  stable  marching  step  size.  ul-cs-P-ij 
Uniform  step  size  was  used  in  this  study,  but  it  is  possible  that 
variable  step  sizes  allowed  Bhu+.ta  and  Lewis  to  march  ^ 

successfully!  however,  no  discussion  of  the  step  size  variation  y 

and  its  relationship  to  stability  appeared  in  their  original 
work. 


Thesis  Advisor t 
Title i 
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Professor  of  Aeronautics  and  Astronautics 
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1.0 


Bhutta  and  Lewis'  have  proposed  a  n*«  and  strikingly 
different  ssthod  of  solving  ths  Parabolized  Navisr-Stokss  <PNS) 
aquations.  This  new  mathod  claims  graatar  accuracy  for  tha  aama 
computing  tima.  In  assanca,  tha  original  idaa  bahind  tha  new 
schama  is  to  solva  tha  normal  aystam  of  non-linear  partial 
diffarantial  aquations  with  tha  aquation  of  stata  lncludad  as  an 
additional  aquation  and  pra'asura  as  an  additional  stata  variabla. 


Tha  antira  mathod  will  ba  raconstructad  and  raaxaminad  in 
datail  hara.  Tha  primary  foaus  will  ba  on  tha  classical 
substitution  for  prassura  in  tha  momantum  and  anargy  aquations, 
and  tha  prassura  rola  in  Rafaranca  1  whara  it  is  kapt  saparata  in 
tha  stata  aquation.  Specifically,  a  convantional  tvo-dimanslonal 
stata  vactor  of  tha  forms 


is  augmantad  in  tha  naw  PNS  schama*  and  appaars  as 


Although  aquations  and  unkowns  consiatantly  incraasa  by 


ona,  two  of  tha  rasulting  stata  alaments  ara  aquivalant  according 
to  tha  stata  aquation.  Tha  aquation  system  is  of  a  mixed 
diffarantial /algebraic  type.  Moreover,  Bhutta  and  Lewis 


explicitly  rule  out  the  need  for  a  sublayer  assumption.  Previous 
PNS  coc'zs  have  introduced  pressure  assumptions  in  the  subsonic 
portion  o£  the  boundary  layer  in  order  to  eliminate  elliptical 
constraints  on  marching  downstream  and  the  development  of 
departure  solutions. 

The  PNS  Equations  were  developed  to  save  on  the  large 
storage  requirements  needed  for  the  full  Navier-Stokes  equations. 
Their  development  arose  from  the  need  to  solve  large  numbers  of 
problems  where  viscous  contributions  are  dominant  in  the 
direction  normal  to  the  streamwise  direction.  In  other  words, 
flows  that  are  boundary-layer-like  are  suitable  to  be  solved  by 
the  PNS  equations.  This  Includes  a  large  class  of  high  speed 
flows  which  are  of  current  interest. 

Previous  PNS  codes  have  used  some  sort  of  sublayer 
pressure  assumption  to  enable  a  marching  solution  in  the 
streamwise  direction.  Lin  and  Rubin* 4  have  used  pressure  from 
experiments  and  have  set  the  subsonic  streamwise  pressure 
derivative  equal  to  the  derivative  at  the  edge  of  the  subsonic 
layer.  Lubnrd  and  Helliwell3  *  1  4  use  a  backward  difference  for 
the  streamwise  pressure  derivative.  Vigneron  et  al3 • 1 4  treated 
the  pressure  derivative  exactly  in  the  supersonic  region  defined 
by 

M«  >  1/  <2->) 

For  other  Mach  numbers  the  pressure  derivative  is  suppressed  by 


the  factor 


( ' 


Cl/<?-l>  3£l-l/(lf(?-l>M« ) 3 

Schiff  and  Stager*  specified  the  subsonic  pressure  to  be  equal  to 
pressure  at  the  bottom  of  the  supersonic  region.  Kaul1 4  uses  a 
global  relaxation  over  the  entire  domain  while  still  maintaining 
the  marching  scheme.  Lighthill11  in  1953  found  that  departure, 
or  separation-like  solutions,  are  observed  in  the  boundary  layer 
equations  whon  the  streamwise  pressure  gradient  is  not  specified. 
Similar  behavior  has  been  investigated  in  the  PHS  equations  by 
Barnett. * 


Figure  1. 1  Supersonic  Velocity  Profile 

Since  the  flow  within  some  region  near  a  surface  is 
always  subsonic  (Fig.  1.1)  and  hence  elliptic,  those  assumptions 
have  been  required  to  change  the  physical  description  from  a 
mixed  hyperbolic/elliptic  to  hyperbolic  character.  Reference  1 
on  the  other  hand,  claims  that  without  changing  the  description 
of  the  physical  nature  of  the  f lowf ield,  i. e. ,  without  sublayer 
assumptions,  the  code  can  be  marched  in  hyperbolic  fashion. 
Justification  rests  there  upon  a  stability  analysis  which  is 
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based  on  a  modified  state  equation.  Further  analysis  will  also 
be  presented  in  the  stability  chapter  of  this  thesis. 

Although  cast  in  the  f ramevork  of  a  PNS  problem,  the 
implications  of  such  modification  of  an  equation  system  to 
achieve  stable  marching  would  impact  in  many  areas  besides  fluid 
dynamics.  Being  able  to  mathematically  change  the  character  of  a 
system  of  equations  without  changing  the  physical  aspects  of 
particular  problems  would  benefit  most  if  not  all  of  the  physical 
sciences  with  promises  of  greater  accuracy  and  efficiency  of 
solutions.  It  is  therefore  of  some  importance  to  rigorously 
examine  the  suggested  new  method. 


2.0 


PRPBLStt  BiSEBIEIIflM 


In  reproducing  the  new  PNS  scheme,  the  emphasis  here  is 
on  the  role  of  the  Individual  procedural  components  of  the 
algorithm:  normalization,  initial  conditions,  step  size, 

iterative  technique,  matrix  formulations,  smoothing,  etc.  and  how 
these  relate  to  convergence  and  marching.  This  study  is 
concerned  more  with  the  evolving  solution  rather  than  the  speed, 
accuracy,  or  efficiency  of  ‘the  scheme. 

In  order  to  understand  how  the  method  works,  the  first 
objective  of  this  thesis  was  to  devise  a  faithful  code  which 
would  provide  sample  calculations  for  a  simple  but  meaningful 
physical  problem.  A  wedge  flow  at  M-3  and  15  and  Reynolds 
numbers  ranging  from  4x10*  to  IxlQ7  are  used.  The  second  purpose 
involves  stability  questions* ••• 1 *’* 4 •* ■  and  includes  eigenvalue 
analyses  of  a  previous*  scheme  and  the  new  PNS  scheme,  as  well  as 
researching  the  behavior  of  systems  of  partial  differential 
equations. 


2. 1  PROBLEM  AS  TESTED 

Reference  1  used  the  new  PNS  scheme  on  a  blunt  body  at 
Mach  25  at  Reynolds  numbers  of  2.92x10s  and  1.72x10*  baaed  on 
nose  radius.  The  reason  for  these  values  was  to  test  the  scheme 
at  high  velocities  and  low  Reynolds  Number  with  large  viscous 
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effects,  where  conventional  PNS  schemes  break  down. 

Since  the  current  study  is  concerned  primarily  with  the 
validity  of  the  concept  in  question,  it  uses  a  simpler  problem t 
a  simple  wedge  flying  at  Mach  3  and  IS  and  Reynolds  numbers  from 
4xlOa  to  lxlO7 .  Reynolds  number  can  be  varied  to  simulate  more  or 
leas  viscous  conditions. 

The  figure  below  shows  the  coordinate  system  of  the  wedge 
and  Table  2. 1  shows  the  freestream  conditions. 


0^  a  wedge  angle 
/3  “  shock  angle 


Figure  2. 1  Wedge  Coordinate  System 
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FREESTREAM  VALUES*  80, OOO  FEET 


DENSITY 
PRESSURE 
TEMPERATURE 
SONIC  SPEED 
VISCOSITY 


8.  6  x  10“  • 
sa.  125 

aas. 99 


977.6 
9.  7  x  10- * 


SL/FT 
LB/ FT* 

•R 

FT/ SEC 
SL/FTSEC 


Tabls  2. 1  Fr»«itr«in  Values 


3.0 


CONSERVATION  EQUATIONS 


As  the  name  suggests,  the  Parabolized  Navier-Stokea 
Equations  are  derived  from  tha  full  Navier-Stokes  Equations.  In 
a  somewhat  similar  approximation  as  for  tha  boundary  layar 
aquations,  only  soma  of  tha  viscous  tarms  ara  ratainad.  Tha  PNS 
aystam  is  applicabla  to  suparsonio  flow  and,  sinca  it  is  valid  in 
both  tha  viscous  and  inviscid  portions  of  tha  flow,  tha 
interaction  between  such  rag Iona  is  includad  automatically.  * 

In  a  normal  suparsonio  viscous  flow,  tha  no-slip 
condition  at  tha  surfaaa  implias  that  tha  boundary  layar  is  at 
soma  point  subsonic  and  transitions  to  suparsonio  flow  within  tha 
boundary  layar  (Figura  1.1). 

A  consaquanca  of  this  subsonic  sublayar  is  that  tha 
atraamwisa  flow  is  alliptic  daapita  tha  outar  suparsonio  flow  and 
tha  solution  cannot  ba  marchad  downstraam  in  a  hyperbolic/ 
parabolic  fashion.  Tha  constraint  introduced  by  sublayer 
assumptions  is  control  of  tha  upstream  communication  of  pressure 
disturbances  through  tha  subsonic  region.  By  making  appropriate 
pressure  assumptions,  such  as  constant  pressure  across  the 
sublayer,  tha  solution  can  ba  marchad  in  tha  straamwise 
direction.  This  has  been  demonstrated  in  a  number  of  previous 
PNS  schemas. *•••*•  14 
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The  novelty  of  the  current  method  la  that  no  auoh 
aaaumptlon  la  wada.  The  olalm  la  that  by  treating  the  equation 
of  atate  aa  a  aeparate  but  coupled  equation,  the  elliptic  nature 
of  the  problem  la  circumvented.  Although  no  aaaumptlona  may  mean 
greater  accuracy.  It  la  neceaaary  to  demonatrate  that  the  atate 
equation  treatment  somehow  altera  the  character  of  the  equation 
ayatem.  Thla  point  will  be  conaldered  further. 


3.1  DERIVATION 


In  two  dimensions,  the  Navier-*Stokes  Equations  in  vector 
form  appear  asx 


q,  *  ♦  CEi  -  Ev],  *  ♦  CGI  -  Gvl,,  -  H  «  0  3.1 

where  q  la  the  atate  vector,  the  comma  indicates  partial 
differentiation  and  E  and  G  are  the  X  and  Z  flux  vectors.  The  1 
and  v  refer  to  inviscid  and  viscous  and  all  variables  are 
dimensional.  Tho  vector  H  contains  the  atate  equation. 

The  X  and  Z  directions  correspond  to  the  components  along 
and  normal  to  the  wedge  surface  of  the  current  problem  (Fig. 

2. 1). 


I 

I 
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Hara  t 


of  atat.a.  Tha  danaity-tsmparatura  product  haa  baan  uaad  hara  for 
conaiatancy  vith  tha  eonvantion  in  Rafaranoa  1. 

Tha  darivatlon  of  PNS  aquatlona  involvaa  an  ordar  of 
aagnituda  analyaia,  auoh  Ilka  tha  boundary  layar  aquatlona.  Tha 
aoat  ooMvon  fora  (Lubard  and  Halliwall,  1973,  1974  >  *  la  obtainad 
by  asauaing  ataady  flow  and  that  tha  atraaawiaa  viaooua 
darivativo  tarva  ara  nagligibla  conparad  to  tha  noraal  (and 
trtnavaraa,  if  3-D)  viaooua  darivativc  taraa.  In  othar  words, 
tha  PNS  aquatlona  ara  darivad  simply  by  dropping  all  vlacous 
taraa  containing  partial  darlvativaa  vith  raapact  to  tha 
atraaaviaa  diraction.  With  thia  in  aind,  it  la  rathar  aiapla  to 
raduca  aquatlona  3.1  to  PNS  forat 

Ei,  «  ♦  Gi,  i  ■  Gv,  i  ♦  H  3. 3 

Tha  atata  vaotor  raiaina  tha  aana  aa  in  aquation  3. 2. 
Sinca  tha  aquation  of  atata  ia  a  aaparata  aquation  and  ia 
algabraic  in  form,  it  ia  tha  aola  contributor  to  tha  aaparata 
vactor,  H. 
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0 


Y  -<*T  I 

The  above  flux  vectors  remain  dimensional. 


Following  the  oonvention  of  Reference  1,  the  non~ 
dleenelonel  variables  becoaet 


■  u* /•’ •  (•  now  indiaates  a  dimensional  quantity) 

*  *’  /a’  • 

•fi'  /P'm 

•  V  /T  m 

«  P'  /  ( Pm  mm*  >’ 

»  p’  /p’  • 

»  X’  /L* 

-  2*  /V 


When  these  definitions  are  substituted  into  the  dimensional 
equations  3. 4,  the  result  is  the  nondimenslonal  PNS  equations: 


Ei,  i  ♦  Gi,  i  ■  €Gv,  *  ♦  H 


The  nondimoalonal  at  at*  aquation  oontalnad  In  tha  H  vaetor  la 
yp  *  ft T.  Unlaaa  atatad  otharwiaa,  all 

aquatlona  and  paraaatara  hanoaforth  will  ba  non-diaansional. 

Tha  factor  €  raaulta  froa  tha  normalisation  prooaaa  and 
la  daflnad  aa 

€  -  Na/Rai  3.6 

Ra«.  la  tha  Raynolda  nuabar  baaad  on  tha  r&faranca  langth,  L*  . 

Rak  **  (aaumL/pa)’  3.7 

For  tha  vlaaoua  wadga  flow  In  thla  atudy  a  oharaatarlatlo 
acala  ^aramatar  la  impliad  by  tha  boundary  layar  d.-*  aplacamant 
thloknaaa,  6* .  In  dimanaional  fora* 

*•’  *  il.720a)X' /<Rat  )*  '•  3.6 

Squaring  both  aldaa,  multiplying  both  aldaa  by  X* /(L*  >*  and 
aiaplifylng  givaa  (in  nondiaanaional  form) 

<&*>*  *  <2. 96)X/RaL 


L*  la  tha  rafarance  langth  unoaan  to  ba  tha  wadga  aurfaca 
location  at  which  X“(X'  /L*  >■!,  ao  that  Ra«.  followa  from: 


R»u  «2.  96/  ( 6*  )• 


3.9 


*•  la  tha  nondlaanaional  boundary  liy*r  diaplacaaant  thlaknaaa  at 
X«A.  Thla  la  aatlaatad  by  knowing  tha  boundary  layar  thlaknaaa 
at  tha  atarting  point.  If  tha  atartlng  plana  la  at  Xal  thara  la 
*  eorraapondlng  Bat  and  L*  for  any  ohoioa  of  1*  and  vloa  varaa. 
for  axaapla,  at  Ha  ■  3  and  with  a  7*  (half  angla)  wadga,  tha 
ahook  angla  la  about  24*.  Proa  Plgura  2.1,  with  X  ■  1, 
Z"tan(ahoek  angla  -  wadga  angla)  or  Z  ■  tan  (24*  -7*)  *  0. 3. 
Aaauaing  a  boundary  layar  that  ooauplaa  10X  of  tha  ahoek  layar, 
ita  thlaknaaa  will  ba  about  0. 03  and  a  raaaonabla  approxiaation 
for  &*  aight  ba  0.01.  Uaing  3‘.  9* 


3.3 


GENERAL  COORDINATE  TRANSFORMATION 


The  PNS  equationa  in  the  form  of  3. 3  are  appropriate  if 
the  grid  i a  everywhere  orthogonal  and  reatangular. 

Unfortunately,  a  aatoh  of  grid  with  the  wedge  and  ahook 
boundariea  requirea  a  tranaforaation.  Nora  generally  the 
equationa  auat  be  prepared  to  allow  a  general  coordinate 
tranaforaation  from  a  nonorthogonal  phyaioal  plane  to  an 
orthogonal  computational  plane.  The  detaila  of  a  apeoifio  grid 
tranaforaation  for  the  wedge  will  be  explained  in  Chapter  4)  the 
general  tranaforaation  of  the  PNS  equationa  will  be  diaouaaed 
here. 


In  eaaenae,  the  goal  ia  to  relate  the  non-orthogonal  X,  Z 
eyatem  to  an  orthogonal  f,  4  eyetea.  We 
take  f  ae  the  etreamwiae  direction,  and  {  ae 

the  croaaflow  direction.  To  begin,  atart  again  with  the  full 
non-dimenaional  ateady  Navier-Stokea  equationa t 

CEi  -  €Ev ] , i  ♦  CGi  -  €Gw], ,  -  H  -  0  3.11 

To  tranaform  to  the  f,  t  ayatcm,  the 

derivativea  (,X)  and  (,Z)  auat  be  changed  to  a  combination  of  the 

derivativea  (,f)  and  <,{>.  By  the  chain 

rule. 
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3.  12 


<,x>  »  <$•,  ,  >  <,  §)  ♦  <4,  * )  <,  4> 

-  in  < ,  *>  *  &  < ,  £) 

<,Z>  »  in  <,i>  ♦  4.  <,4> 

This  assumes  a  most,  general  transformation  of  the  form: 

i  •  ax,z) 

4  *  4<X,  Z) 

Now  use  3. 12  in  3. 11  to  get 

C!x  (Ei-€Gi>  it  (Ev-€Gv>  3,  i 

♦  C4,CEi-€Gi}  4i  <Ev-€Gv>  3, 4  *  0  3.13 

In  keeping  with  the  assumptions  to  derive  the  PNS 
equations  from  the  Navier-Stokes  Equations,  the  streamwise 
viscous  derivatives  are  omitted,  which  leaves: 

C in  (Ei-€Gi>  ],f 

♦  C4x  (Ei-€Gi>  ♦  4i  {Ev-€Gv>  1,4  3.12a 

The  Ev  and  Gv  viscous  vectors  contain  X  and  Z  derivatives  which 
must  also  be  transformed  to  £  and  4>  The 

entire  derivation  may  be  found  in  Appendix  A.  The  transformed 
PNS  equations  become: 

F,  ,i  *  F,  ,  4  ■  €S,  4  ♦  H  3.  14 

Here  the  notation  is  different  to  emphasize  that  the 
transformation  has  been  completed  and  H  is  still  the  vector 
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containing  the  aquatic n  of  atata  terras.  ?he  H  vactor  ia  not 
tranaformed  bacauaa  of  ita  aigabraic  nature.  Qthar  vactors  are 


S 


And: 


/ 


(1/J) 


\ 


pU, 

puUt  *  ft  P 
PwUi  ♦  it  p 
( T/ (7-I )  ♦V*  /2>pU, 
0 


(1/J) 


put 

PuU«  ♦  4>  p 

P*U«  ♦  &  p 

<T/<7-1)+V«/2)pU, 

0 


(p/J) 


M„u,  4  ♦  (Hi«u,  4+M*  i  wf  4)  /3 

M»  *»  S  ♦  ( M*  t  u#  {+Si  i  {)  /3 

Mo  ( T,  4/  ( Pr  <  7- 1  >  >  +uu,  &+**,  4) 

♦  <M*  i  uu,  4*M«  (  »¥,  4 

*M*  i  <  wu,  4+uw,  4)  >  /3 


Ui  ■  ft  u  ♦  4*  w 

Ut  ■  4k  u  ♦  4>  w 

fit,  ■  4»“ 

Mi  i  =4*4* 

M, ,  -  4ia 


3.  IS 
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Mo  ■  M*  *  +  M*  i 
V*  »  u*  ♦  w. 

J  *  Transformation  Jacobian  defined  in  Eq.  (4.7) 
The  Jacobian  is  included  to  make  the  entire  transformation 
conservative. s  This  will  be  seen  more  clearly  in  Chapter  4. 


Equations  3.  14  are  to  be  solved.  However,  in  preparation 
for  a  later  stability  analysis,  a  differential  version  of  the 
equation  of  state  introduced  in  Reference  1  must  be  noted: 

8<P,  §  *  P,  5)  ♦  vp  *  pT 

Here  8  is  a  "small"  number.  With  this  version  of  the  equation  of 
state  the  final  versions  of  Ft  and  F«  become: 


F, 


F. 


<1/J> 


<1/J> 


, 


Pi J. 

<ouUt  ♦ftp 
(OwUt  ♦  ft  p 
(T/('>“1)+V*/2>^U 

0p 


/“  \ 

pull,  ♦  ftp 

i“wU«  ♦  ft  p 

<T/<7-l)+V«/2)(aU. 

\  ®P  / 


3.  16 
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4.  O 


GRID  AND  GEOMETRY  ANALYSIS 

The  viscous  wedge  problem  clearly  involves  a  boundary 
layer.  The  presence  of  viscosity  implies  that  the  location  of 
the  shock  is  unknown  beforehand,  unlike  the  inviscid • problem. 

The  shock  locus  is  found  as  the  solution  is  marched  along  the 
wedge.  For  a  shock  fitting*  approach,  the  shock  serves  as  the 
upper  boundary  of  the  described  domain.  A  convenient  grid  is  one 
that  conforms  to  the  physical  boundaries. 

In  order  to  solve  the  flow  field  numerically  consider  a 
rectangular  grid  in  a  computational  plane.  The  relationship,  or 
transformation,  between  the  physical  plane  and  the  computational 
plane  provides  the  metrics  that  were  developed  in  Chapter  3. 0, 
equations  3.12: 

4*  ,  ,  5*  ,  & 

To  perform  the  metrics  analysis,  begin  with  a  general 
relation  between  the  physical  and  computational  planes. 

This  general  relation  has  the  form: 

i  *  *<X, Z>  4.  1 

$  -  ax,  2)  4.2 
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and  serves  as  a  mapping  between  one  plane  and  the  other 


From  the  chain  rula 


( ,  f )  *  <  X,  fX  ,  x )  -X  Z,  f X  ,  z )  4.3 

<,£>»<X,£X,x>XZ,£X,z>  4.4 

Sine*  tha  transformation  ia  from  tha  physical  to  tha 
computational  planes,  tha  derivatives  (,x>  and  (,z>  must  be 
written  in  terms  of  the  derivatives  ( , f )  and 
(,£).  Solving  4.3  and  4.4  simultaneously  gives: 

< ,  x )  "  JC  <  Z,  £)  (,£)-(  Z,  f)  < ,  £>  3  4.3 

(,  z)  *JC  ( X,  f)  (,£)■*  (X,  £)(,£)  1  4.6 

where  the  Jacobian  of  the  transformation,  J,  is  represented  by 

J»1/C(X,  *XZ,£>-<X,  £XZ,*>3  4.7 

and  appears  in  equations  3.  15  and  3. 16. 

Now  the  general  transformation  of  Chapter  3  can  be 


completed.  By  comparing  4.  5  and  4.  6  to  equations  3. 12  it  is  sei 
that  the  required  metrics  in  terms  of  the  physical  plane  are: 


f,  -  “X,  {( J> 


4.9 


&  •  -Z,  f<J)  4.10 

Cr  «  X,  f<  J>  4.  11 

Sine*  ths  physical  gsomstry  is  nonuniform  vhils  ths 
computational  plans  rsmains  constant,  thsss  mstrics  must  bs 
rscalculatsd  at  sach  vsrtical  <L)  nods  at  all  strsamwiss  <J> 
locations.  Sss  Figurss  4. 1  and  4. 2. 

Ths  final  contribution  to  ths  grid  transformation  rslatss 
to  fixirg  A£  and 
AC  sines  thsy  ars  ussd  to  find 

dsrivativss  in  ths  computational  domain.  At  sach  strsamwiss 
station  thsrs  ars  always  ths  sams  nuvnbsr  of  points  vsrtically, 
say  LNAX.  Sines  S  goss  from  0  to  1,  (Fig.  4.2). 

AC  «  1/ ( LMAX-1 ) 

Ths  rangs  of  Z-coordinatss  in  ths  physical  plans  vary  with  shock 
laysr  thicknsss,  but  C  is  always  mappsd  to  ths  rsgion  0 
to  1  with  squally  apacsd  points  AC 
apart. 
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The  f  coordinate  goes  from  lb  to 


fc « « .  dependent  on  the  atart  and  end  of  the  physical 
problem  and  the  transformation.  Therefore! 

U  Xfbo.-lb  > / < JMAX-1  > 

Where  JHAX  is  the  number  of  streamwise  stations,  the  initial  data 
station  being  station  number  1. 

a 

Ah  more  node  points  are  added  vertically  and 
horixontally.  the  accuracy  improves  according  to  the  accuracy  of 
the  governing  finite  difference  equations  (Chapter  S).  This  is 
in  contrast  to  Reference  1  in  which  the  computational  grid 
appears  to  hold  Af  and 
At  constant  at  1.0. 

It  remains  to  calculate  the  metrics  using  4. 8  -  4.11. 

One  could  do  so  numerically  by  taking  differences.  For  example. 
(Equ.  4.10; 


&  »  -AxJ/Af 

*  “ J (Zj.i  ,  l  "Zj  ,  L 

Where  J  muat  also  be  found  by  differencing.  Of  course  knowing 
the  physical  grid  locations  at  the  j  ♦  1  station  in  order  to  do 
the  differencing  requires  knowledge  of  the  shock  location.  The 
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I 

j 


\ 

I 


■hook  prediction  will  be  dieoueeed  in  the  chapter  on  boundary 
condition*. 

If  one  knowe  the  epecific  correspondence  between  the 
physical  and  computational  planes,  ae  in  equations  4. 1  and  4. 2, 
the  metrics  may  be  found  explicitly.  For  the  wedge  problem,  the 
grid  transformation  is  given  by* 

f  -  X  4.  12 

{  *(l/s)sinh~  * (sinh(s) CZ/Z«mocn 3  >  4.  13 

Equation  4. 13  was  chosen  to  cluster  the  grid  points  in  the 
boundary  layer4.  With  an  "a1*  value  of  3.0  for  example, 
approximately  40X  of  the  grid  points  lie  within  the  boundary 
layer.  The  "s"  value  can  be  varied  to  give  the  desired  grid 
clustering.  Equation  4.13  can  be  rewritten  as: 

Z  *  <sinh(s()  /ainh(o)  >  CZ«N  (x)  3  4.14 

Where  Z« M  is  a  function  of  X  and  hence  (by  4. 12),  of  f 
also.  By  4.12,  4.13  and  4.  14* 


X.  t  -  1 

x,«  -  0 

Z,  t  •  (einh(s{) /einh(s)  >  CZ’  •*  3 


Z,  $  ■  (■ )  <  Zt  m  >  ( cosh  ( s£)  /ainh  ( a  >  ) 


vhara  subscript  SH  rsfsrs  to  ths  vslus  *t  ths  shock  snd 


Z*  ■  M  ■  (  Zj  «  I  ~Zt  )  ■  N  /  At 

Pros  ths  sbovs  formulas  and  aquations  4.8  -  4. Ilf 
fi »  ft  i  {■  i  and  fr 

can  bs  constructad  at  aach  straasvisa  station  and  at  aach 
’-artical  r. oda  as  tha  solution  advancas  in  tha  straamvisa 


direction 


km  wa  hav*  aaen,  tha  non-dimanaional  2-D  PNS  aquationa 
can  b*  writtan  mat 

Pi  *t  *  Fa#  *  *  €Sf  «  *  M  S.  1 

whara  F» ,  F. ,  S,  and  H  ara  vactora. 

Tha  nav  achama  la  implicit  and  itarativa  which  maana  that, 
tha  aolutlon  la  aovad  forward  in  apaca  from  indax  j,  aa y,  to  j+1 
and  than  ltaratad  in  "paaudo  time"1  from  n  to  n+1  until  it 
convargaa  (Fig  S. 1 ) .  Tha  tarm  "paaudo  tima*  rafara  to  tha  fact 
that  only  tha  oonvargad  ataady  atata  aolutlon  haa  phyaical 
meaning. 


n*time  indax  L«vartical  indax  j "horizontal  indax 


Figura  S. 1  Index  Notation 


With  referenee  to  Figure  9. 1,  the  converged  solution  st 
stops  j  snd  j -I  srs  uasd  ss  ths  starting  solution  for  stop  $*1. 

At  stop  j  ♦  1,  ths  solution  is  thsn  itsrstsd  until  oonvsrgsnos.  At 
stop  j+l  vs  ha vs i 

<F,,f  ♦  F.,«  ■  S,S  ♦  «>■.**»••♦»  5.2 

Sinos  ths  msthod  is  implicit,  all  vsrtioal  (.-nods  grid  values  at 
ths  n+1  Isvsl  ars  found  simultaneously  in  oontrast  to  an  sxplioit 
scheme  vhich  would  find  ths  n*l  values  one  at  a  time  from  ths 
already  known  n  values. 

Ths  above  equation  is  not  useful  since  the  solution  at 
n+1  is  unknown  --  the  n+1  level  must  be  tied  to  the  n  level.  If 
we  assume  that  the  solution  at  level  n*l  is  close  to  that  of  the 
n-th  level1 ,  Taylor  Series  expansions  in  pseudo  time  for  each 
term  in  9.2  gives i 


p,  j  t  «  F,  J**.«  ♦  <A,  **  1  •  »  1  )  5.3 


F,J**.»*‘  ■  F«  ‘  "  ♦  (A.  )  5.4 


SJ  ♦».<*♦»  -  SJ  ♦  ‘  "  ♦  <H**»  •  >  5.5 


HJ  ♦  <A0J**' )  5.6 
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Where i  £q" • 1  ■  qj  ♦ « •  »  * »  -  qJ ♦ ‘  »  and  theae 
expanaiona  ara  valid  at  each  vartloal  node. 


An  A> ,  Aa »  and  H  ara  callad  tha  Jacobian  matricaa,  not  | 

to  ba  oonfuaad  with  tha  tranaformation  Jacobian  (Chapter  4). 

Thay  ara  S  x  S  matricaa  ainca  thara  ara  five  elementa  in  tha 

atata  vactor.  Tha  Jacobian  matricaa  ara  formed  by  taking  tha 

partial  derivativea  of  tha  flux  vactora  vith  raapact  to  tha  atata  j 

•  i 

vactor.  Tha  alamanta  of  tha  Jacobian  matricaa  ara  givan  in  i 

Appandix  B.  I 

v 

Subatituting  aquationa  S.  3  -  5.6  into  5.2  gives*  | 


CF,  •£q***  3,f 

♦  ♦  A.J*‘«  •  Aq»*‘ 3,4 

■€ESJ**»"  ♦  MJ  ♦  *  .  »  Aq"  ♦  t  1,  { 

♦  CHJ  *  ‘  •  "  +  A*  •>  *  1  •  "Aq"*‘  3 


5.7 


Tha  atraamviaa  derivative  was  givan  spacial  traatmant  in 
Rafaranca  1  ao  that  tha  truncation  error  vould  ba 
0( £q" *  *  )• 

inataad  of 

0<£qJ  ♦  ‘  )■ 

aa  in  conventional  PNS  achamaa. •  •  •  •  1  *  •  *  4  •  ‘ *  It  ia  preferable  to 
have  temporal  arrora  ainca  they  vaniah  for  tha  converged 
aolution. 
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The  streamvise  derivative  from  5.2  becomes: 


/ 


(F,  ,?>**  1  •  ■*  1  »  (FtJ*1-"*1  -  0<A?> 

«tA,  *  *  ‘  •  •  Aq*  *  »  +  (F»  ‘  "  -Ft  J  >  3/A?  ♦  0<A?,  (  Aq*  *  1  >•  > 

■At  ■»  ♦  ‘  •  "  Aq"  ♦  1  /  A?  ♦  <  F»  ,  ?>*  ♦  *  •  " 

With  this  formulation  of  the  streamvise  derivative,  equation  5. 7 
appears  as: 

C  <A,  /  A?  -  Aq  )  (  A  ) 

♦  <  At  -€M ) ,  £(  Aq"  *  *  )  3 1.  J  *  *  •  " 

■  ~ C Ft ,  ?  +  F.,«-  €S,£-  H3L  J  ‘  •  " 

■  Gu  J  ♦  ‘  "  ,  say  5.  8 

Equations  5. 8  and  5. 7  are  equivalent  but  5. 7  has  been  changed  to 
explicitly  show  how  the  streamwise  differencing  is  performed. 
Equation  5.8  is  accurate  to  □(Aq"**>>  in  pseudo 

time  and  is  conservative  in  the  limit  of  convergence. 

It  is  convenient  to  call  the  right  hand  side  of  equation 
5.8  0'.'**’"  for  ease  of  notation.  It  is  seen  that  on  the  right 
aide  all  values  are  in  terms  of  the  n-th  level,  which  are  knovn. 
On  the  left,  the  coefficients  of  Aq"**  are  also 
from  the  n-th  level.  The  problem  is  elliptic  in  the  £ 
direction  so  central  differences  are  used  for  vertical 
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differencing.  Upon  using  a  central  difference  on  the  left  of 
5.fi,  the  final  form  becomes: 


<  At  /At  -  Ao  > u  **  1  •  "  <  Aq"  *  ‘  ) 

» 

♦  {  <  (A.  -€«>,.♦  , /(2A5)  )  <  Aq"*  ‘  ) 

-  <<A.-€M)l_i /(2A£>  HAq**1  >>•»**•  " 
■Gl  J  *  1  •  * 


5.  9 


Looking  at  5.  9,  the  block  tridiagonal  form  begins  to 
appear.  On  the  right,  there  is  the  vector  G-*  ♦  1  •  "  at  all 
vertical  L  nodes.  On  the  left,  there  are  three  matrices;  one  at 
node  L,  one  at  node  L+i  and  the  other  at  node  L-l.  The  equation 
solves  for  the  vectors  AqL «* ♦ ‘ .  Once 
Aqt."*1  has  been  found  with  a  block  tridiagonal 
solver,  the  new  state  vector  can  be  found  from: 

qtj  t  *  ♦  Aqt.  j  ♦»,»*♦  »  5.10 

Equations  5. 9  and  5.  10  do  not  address  smoothing  considerations, 
which  will  be  added  in  Chapter  7. 


To  illustrate  the  block  trldiagonal  construct  a  bit  more 


clearly,  aet 


A  *  “C<A«-€M)/(2A6)3l.-i  J*‘*  " 


C  »  t  ( A«  -6M )  /  <  2A£>  ]u  *  »  J  *  1  •  "  5.11 

B  *  (A,  /L§  -  A0  »*  1  •  " 


Then  the  block  tridiagonal  equation  appears  as: 


BC 

ABC 


I 


\ 


X 


/ 


Gl  «  a 


\ 


5.  12 


ABC 
AB 

Boundary  conditions  for  the  matrix  operations  will  be  explained 
in  the  next  Chapter. 

A  remaining  question  is  how  to  find  the  right  side  of 


4qiL-i 


Gt,  ■ 


uni-1 


3.8;  namely,  the  derivatives: 


As  mentioned  earlier,  the  £  derivatives  are  modeled  by 
central  differencing,  i.e. , 

Ft ,  £  -  "/<2A£> 

S,  £  ■  (  Sl  ♦  »  “Sc  -  t  >**♦*•  "/<2A£) 

The  streamwise  derivative  is  a  one-point  backward 
difference  as  seen  in  the  d'eveiopment  between  S.  7  and  5. 8: 

Ft  *  t  3  (Ft-1*1*  --F»  J  ),./££ 

In  all  such  differencing,  the  metrics  are  also  included.  In 
differencing  the  viscous  terms  the  viscosity  coefficient  is  also 
differenced.  Viscosity  is  determined  by  Sutherland's  Law  in 
dimensional  form:9 

p  »  B(T9/« ) / ( T  ♦  S)  5. 13 

where  T  is  dimensional  temperature  and 

B  ■  7.3025  X  10- T  Ibf  /  <  f  tsecR*  7  •  ) 

S  »  198. 72° R  for  air 

The  Llock  tridiagonal  form  of  equation  5. 12  is  solved 
with  a  subroutine  that  was  originally  written*  for  an  academic 
subject  requiring  the  solution  of  a  block  tridiagonal  matrix  made 


up  of  3x3  matrices  beginning  at  node  1.  It  was  modified  here  to 
allow  5x5  matrices  and  to  start  at  node  2. 

The  block  tridiagonal  form  (5.12)  in  Reference  1  was 
solved  using  stored  forms  of  the  inverse  matrix.  The  scheme  of 
Reference  1  computes  new  matrices  in  the  block  tridiagonal  form 
only  at  the  first  iteration  of  each  streamwise  step.  The  current 
code  computes  these  matrices  anew  for  each  iteration,  before 
solving  the  block  tridiagonal  matrix.  This  revised  procedure  does 
not  affect  the  converged  solution  but  may  change  the  number  of 
iterations  required  to  obtain  the  converged  solution. 

The  tridiagonal  solution  advances  the  solution  in  pseudo 
time  only,  from  level  n  to  level  n+1.  Once  convergence  is 
achieved,  the  solution  marches  from  level  j  to  j+1  by  some  method 
of  prediction.  Far  example,  either  the  newly  converged  j 
solution  becomes  the  first  j+1  solution,  or  using  solutions  at  j 
and  j-1,  extrapolation  predicts  the  first  solution  at  j+1.  The 
iterations  take  place  at  one  streamwise  location  rather  than 
globally  aver  the  whole  field  as  in  previous  iterative  PNS 
schemes.3’  1  *  This,  along  with  the  separate  equation  of  state, 
makes  the  current  scheme  significantly  different  from  past 
algorithms. 
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6.0 


BOUNDARY  CONDITIONS 


There  are  three  distinct  boundaries  for  the  wedge 
problem:  The  wedge  surface,  the  shock  wave,  and  the  initial  data 

plane.  Each  will  be  covered  separately. 

L 

I  6. 1  SURFACE 

f 

I 

Five  boundary  conditions  are  required  at  the  wall 

i 

f 

corresponding  to  the  five  state  variables.  These  conditions 
are:  * 

i 

1)  Density  consistent  with  equation  of  state 

2)  No  slip  condition  for  u  velocity:  <pu)»  »  0 

3)  No  slip  condition  for  w  velocity:  (<ow>»  ■  0 

4)  Specified  wall  temperature,  T»  >1.0  for  Mach  3 
and  T(  ■  3. 0  for  Mach  IS  flew. 

5)  Zero  pressure  derivative  in  the  body  normal 
direction,  <P,  £)  =  0 

Here  subscript  1  refers  to  a  wall  value  and  condition  #5  is 

Justified  by  a  boundary  layer  type  analysis  performed  at  the 

wall. ‘  The  true  nondimensional  equation  of  state  is 
7P~pT  ■  0.  However,  as  will  be  discussed 

in  Chapter  8,  a  modified  equation  of  state  is  taken  to  be: 

7p  ~  pT  ♦  0CP,  f  ♦  P,  63  ■  0  6.1 
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Because  of  condition  #5,  at  the  surface  this  reduces  to 

<  yp  -  P’C  ♦  0P,  t )  i  -  0  6.2 

Where  8  is  an  "arbitrary”  (small)  parameter.  For  sufficiently 
small  8  the  model  approximates  the  true  equation  of  state  if 
<P,  f)*0(l). 

Equation  S. 9  applied  at  the  first  point  away  from  the 
wall  (L*2)  results  in  three  matrices  on  the  left  and  one  vector 
on  the  right-hand  side: 

A,  Aqt  +  B.  Aq.  *  Ca  Aq,  ■  C*  6.  3 

This  vector  equation  represents  5  equations  for  the  5  state 
variables.  Since  the  subscripts  refer  to  specific  nodes, 

At  Aqt  is  known  from  the  boundary  conditions.  It 

is  not  part  of  the  block  tridiagonal  matrix  which  has  the  form 

(see  Equ.  5. 12) : 
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where:  B*t  *  B«  ♦  f 'A*  ) 


G'uhu-i  *  Q 


l  M  I  -  I 


“  Cl.  max  Aqt 


Sine*  Aqi  *  (q"*‘  -  q"  >,  is  unknown 
(specifically  q**‘),  this  boundary  condition  can  b*  treated 
implicitly,  transformed  to  the  delta  form,  and  combined  with  the 
B«  Aq*  term. 


A  model  for  condition  #5  is  that  the  pressures  are  equal 
at  the  wall  and  the  adjacent  node. 


Pt  -  P. 


It  fallows  therefore  that 


Ap»  *Aps 


where  Ap»  is  the  fifth  component  of 

Aq.  The  elocity  terms  are  also  straight  forward 

since  from  conditions  2  and  3: 


A(  PU  )  i  a  A(pw)t  3  0 


For  the  p T  stat  variable  we  have  (from  S.  2  and  6.5) 


A(pT)t  »  (-ypt  ♦8p»  ,  i)"+ ‘  -  (<ypt  ♦0p»  ,  $>" 

*  '  *-0p«,f)"*‘  -  ( 7P«  ■*’0p*  ,  f)' 

■  A( pT )  • 


j\ 

:>«l 


I 

■<‘l 


And  finally. 


A 0,  »  (^aT)t/Ti  ■  ( AoT ) ,  /  T, 


Combining  aquation*  6. S  -  6.6  gives t 


A(pu)t 


A(pw>t 


fA<pT>,/T, 


A(pT), 


A(  (OT )  * 


The  elements  of  the  At  matrix  which  multiply  these 
components  of  Aqi  can  be  combined  with  the 
elements  of  B«  which  multiply  like  components  of 
Aq« ,  leading  to  B*a  of  equation  £.4. 


The  A«  matrix  (Appendix  B)  is  filled  based  on  boundary 
conditions  1-5.  The  velocity  components  vanish,  the  temperature 
is  given,  and  the  pressure  at  the  wall  is  taken  from  the  node 
directly  above  the  wall.  Density  is  then  found  from  the  modified 
state  equation. 


-  <->pi  ♦  C0<P,  *  *  ‘  •  "  -P,  *  )/Aff3  >/T, 


6.  10 


i 

e-.i  >ru  s-v.-  tt-j 


f 


There  ere  two  parte  to  the  shock  boundary  condition:  The 


location  of  the  shock  relative  to  the  surface  of  the  wedge,  ZaM 


and  the  angle  the  shook  surface  makes  with  the  freestream  flow. 


The  shock  angle  determines  the  state  variable  values 


across  the  shock  and  the  location  of  the  shock  is  used  to 


equalize  mass  flow.  Since  the  flow  is  viscous,  the  boundary 


layer  will  force  the  shock  surface  outward  and  at  the  same  time 


the  shock  angle  may  change  as  a  function  of  streamwise  location. 


The  algorithm  must  allow  for  changes  in  both  angle  and  location. 


The  values  of  the  state  variables  at  the  shock  must  be 


consistent  with  those  in  the  interior  of  the  shock  layer.  In 


other  words,  the  shock  values  must  be  coupled  with  the  interior 


values.  To  preserve  this  coupled  nature,  the  shock  angle  must  be 


obtained  from  points  Interior  to  the  shock. 


One  way  to  do  this  is  to  solve  for  pressux'e  or  density  at 


the  shock  based  on  interior  points  and  then  use  the  Rankine- 


Hugoniot  conditions  to  find  the  shock  angle.  Using  backward 


differencing,  any  value  at  the  shock  can  be  obtained  from  points 


Inward  of  the  shock: 


ft.  nai  3  K  (  f  i_  h  a  «  -  4  ~  2f  l  a  a  i  -  i  “  ft.  a  a  i  -  a  ♦  4ft.  a  a  i  -  l  )  6.  1 1 


where  LltAX  is  th*  nod*  at  th*  shook. 

Us*  6.11  to  find  p«M  or  Am  and  tha  Rankine-Hugoniot 
conditions  to  find  ths  shock  angle,  fl.  This  shook  angls  will 
thsrsfors  bs  consistent  with  stats  variable  profiles  calculated 
at  each  iteration.  Once  P«*  or  Am  are  found 

remaining  variables  are  found  from  th*  Rankine-Hugoniot  equations 
and  th*  state  equation.  Th*  velocities  across  the  shock,  uaN  and 
w« M ,  are  found  from  geometric  considerations  detailed  in 
Appendix  C.  On  the  first  iteration,  before  any  profiles  have 
been  calculated,  either  th*  shock  pressure  or  density  can  be 
extrapolated  from  upstream  values  with  a  simple  Eulerian 
integration.  The  shock  values  found  at  each  iteration  become  th* 
values  at  the  n-»l  level. 


The  global  mass  conservation  procedure  considers  a 
problem  separate  from  the  shock  value  calculation  and  independent 
of  th*  coupled  nature  of  th*  shock  angle.  The  shock  location 
determines  the  maximum  Z  value,  which  determines  the  grid 
distribution,  which  determines  the  metrics.  To  begin  this 
routine,  the  shock  location  must  be  predicted  when  stepping  from 
j  to  j+1  prior  to  any  iterations.  At  each  iteration  th*  ZSM 
position  is  adjusted  to  equalize  mass  flow  vith  freestream  mass 
flow. 

The  method  used  to  predict  the  shock  location  is  based  on 
a  paper  by  Chausee,  et  al. 7  His  method  is  written  for  a  general 
three  dimensional  system  but  here  reduces  to: 
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Z«MJ*‘  *  Z«MJ  ♦  AxttanfiJ > 


6.  IS 


vhara  Ax  la  tha  atraamvlaa  atap  aiza.  Thla 

pradlotlon  la  uaad  only  at  n-1.  At  furthar  itarationa,  Z«M 

adjuatmanta  ara  baaad  on  ovarall  maaa  flow  conaidarationa. 

Tha  aaaa  flow  aquallzatlon  coneapt  vaa  not  part  of 
Chauaaa’a  achaaa.  It  vaa  auggaatad  In  Rafaranoa  1  alnoa 
Chauaaa'a  mathod  raaulta  In  aaaa  flow  arrora  of  t  2. OX. 1  Bhutta 
and  Lavla  auggaat  aovlng  Z% M  until  aaaa  flow  arrora  ara  laaa  than 
±  0. IX.  Datalla  of  both  tha  ahock  pradlotlon  and  tha  maaa  flow 
calculation  ara  praaantad  In  Appandlx  C. 

Raf arrlng  back  to  aquation  6. 4,  tha  ahock  boundary 
concarna  tha  Cl  N«  *  M  *  *  *  *  ‘  vactor.  And  aa  In 

aactlon  6.1,  tha  Clh««  matrix  la  f Iliad  ualng  tha  ahock  valuaa. 
Unllka  tha  wall  boundary  aondltlon,  £qL mi'*1  la 
known.  Tha  nawly  calculatad  ahock  valuaa  at  aach  ltaration 
bacoma  q" * ‘  and  tha  laat  computad  ahock  valuaa  ara  q" .  In  thla 
caaa  AqL»»*  can  ba  formad  explicitly* 

-6qui»«***‘  *  <q**‘  -  q*)i,«i  6.16 

Tha  right  hand  alda  of  6. 4  la  than  corractad  for  thla  boundary 
condition  with: 

0*  l  m  i  •  i  *  Gl  m  •  *  -  »  ~  Cl  mi  (  AcJl  mi**‘  )  6.  17 


SO 


The  Initial  data  plana  (IDP)  ia  tha  starting  point  for 
tha  numerical  solution.  Using  tha  currant  nondimanaionalixation, 
tha  IDP  is  at  X»1.0.  Initial  valuaa  for  all  fiva  stata  vtctor 
componanta  ara  pradictad  or  assumed  from  tha  wall  to  tha  shock. 
Initial  data  which  togathar  with  tha  naxt  initial  pradiction 

battar  satisfy  tha  govarning  aquations  (Equ.  3.8)  rasult  in 
fastar  convarganca  for  tha  first  j+1  location. 

In  Rafaranca  1  a  blunt  body  starting  coda  was  used  to 
obtain  an  initial  condition  for  a  blunt  cona  configuration)  a 
starting  solution  for  tha  prasant  viscous  wadga  flow  was  not 
availabla.  Initial  conditions  wars  based  instaad  on  inviscid 
flow  with  rafinaman  a  basad  on  tha  addition  of  a  boundary  layar 
and  a  mass  flow  calculation. 

Bacausa  of  this,  tha  prasant  initial  conditions  ara 
approximations  at  bast  for  a  "corract"  wadga  flow.  However,  tha 
approximations  ara  baliavad  to  ba  raasonabla  basad  on  an  assumad 
boundary  layar  that  was  10X  cf  tha  initial  local  shock  layar 
thicknass  with  a  spacifiad  distribution  from  tha  wall  to  the  adga 
of  tha  boundary  layar.  Sinca  tha  mass  flow  must  balance,  tha 
shock  location  and  tha  stata  variabla  profilas  ware  than  adjusted 
to  match  tha  fraastraam  mass  flow  with  that  of  tha  IDP. 


Am  m  n»t«urt  of  tha  quality  of  the  initial  data,  the 
right  hand  aid*  of  Equation  3. a  was  monitorad.  Tha  right  hand 
aida  ia  aimply  tha  govarning  PHS  ayatam  tarma  writtan  at  tima 
laval  n.  A  oonvargad  aolution  oorraaponda  to  a  vaniahing  right 
hand  aida  to  aoma  aooaptabla  laval.  Tha  roagnituda  of  tha  largaat 
componant  of  tha  right-hand  aida  for  tha  firat  itaration  aftar 
tha  IOP  indicataa  how  vail  a  conaiatant  aolution  vaa  achiavad 
with  tha  initial  data  and  tha  pradiction  at  j+1.  Aftar  aatting 
up  tha  initial  data,  profilaa  vara  amoothad  to  bland  tha  ragion 
at  tha  top  of  tha  boundary  layar.  Figuraa  6. 1  and  6. 2  ahow 
typical  initial  data  for  Mach  2  and  for  Mach  IS  flow.  Tha  figuraa 
ara  only  a  rapraaantation  ainca  Raynolda  numbar  vaa  variad  and 
tharafore  diffarant  initial  profilaa  vara  uaad  for  diffarant 


caaaa 


7. 0  SMOOTHING 

The  final  form  of  the  finite  different  equations  are 
given  by  equation  5.9  with  the  right-hand  side  as  in  5.8.  For 
simplicity,  this  can  be  written  as: 

CT ]  "  Aqt.  "  *  1  *  CRHS3"  7.1 

where  CT3"  is  the  block  tridiagonal  matrix  and  CRHS]"  is  the 
right-hand  side  vector.  The  equation  as  it  stands  solves  for 
Aq"**  and  as  yet  there  is  no  smoothing. 

As  stated  in  Reference  1  and  by  Schiff  and  Steger, * 
central  differencing  produces  oscillatory  behavior  which  must  be 
damped.  Since  Equation  6.  1  is  second  order  accura-ce  in  the 
£  direction  due  to  the  central  differencing,  it  is 
possible  to  add  a  term  of  0<A£>*  as  a 

smoothing  parameter  without  formally  affecting  the  second  order 
accuracy. * 

Reference  1  chooses  the  form  of  this  smoothing  parameter 

to  be: 


wtf  <qJ  ♦  1  >  3  <Aq  )* 


7.  2 


where  f  is  an  appropriate  vector  and  w  is  some 


constant 


Bhutta  and  Lewis1  chose  the  vector  f  to  bet 

f  *  *EA,q,«  «•  <  A«  q,  fi£) ,  £  -  €<Mq,  $6> ,  «  -  A0  q,  £3  7.3 

and  w  ■  0.  or  1.0  for  no  smoothing  or  smoothing. 

If  7.2  is  added  to  S. 1  we  obtain: 

F, ,  §  ♦  F«  ,  £  •  €S,  £  H 

♦  <*f  <q *  *  ‘  )  <A£>«  7.  4 

Substituting  for  f<qJ*M  (and  using  the  disc. salon  of  Chapter  5) 
results  in: 

EF,  ♦  A»  (-uq,  £(A£)«  >/43->  '  ‘ 

♦  CF«  ♦  A«  ( -«q,  £(  A5) •  ) /4  J  J  *  1 , 5 
*€ES  ♦  M(-«q,  £(  A£)«  J/41-*  '  1  ,  £ 

+  EH  Aq  <-«q,  ££(A£;«  )/43«*  » 

*F,  *  *  <££)*  ) 

7.3 

Now,  define  the  quantity1 

-«q,  6S<  ^  /4  ■  QJ*‘  -  qJ*‘  ■  0<A£)‘ 

7.6 


mo  that  (□■<♦»  -  .  Q<£5>4  7.7 

or,  to  second  order  accuracy: 
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qj*»  «  QJ*»  ♦  «Q,  S6<A5>* /4 


7.  8 


Consider  s  Taylor  Series  expansion  of  Ft  (Q)  around  q: 

Ft  <Q)  *  Ft  <q)  ♦  Ft  ,  q  <0-q>  ♦  Ft  ,  qq  <Q-q)*/2  +  .  .  . 

«  Ft  (q)  ♦  A.  <-<q,  $S<  A£>«  > 

♦  0  </&>♦ 

Similar  constructs  can  be  obtained  for  F* ,  S  and  H  and  in  terms 

of  an  intermediate  solution  Q  and  Equation  7.5  can  be  written  as 

F»  ♦  F.  <Q)4**,«  -  €CS(Q)-»*<  J,  £  *  H(QJ  *  1  )  +  F,  J 

7.9 


The  intermediate  solution,  is  related  to  q  by  Equation  7.8. 
Comparing  Equation  7,9  to  the  development  in  Chapter  5,  we  can 
see  that  7.9  in  terms  of  Q  can  be  handled  like  5.8  in  terms  of  q 
In  other  words,  the  governing  equation  can  now  be  written  as 


CAAfl**1  ♦  B&2"**  ♦  C^Q"*‘3J**  •»  Gt J  *  *  •  " 


where  A, B, C  are  the  tridiagonal  matrices,  or 


7.  10 


CT]"£G»*‘  «  CRHS3" 


7.  11 
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This  is  ths  equation  to  bs  coded  end  solved.  Once 
££3-*  *  *  •  *  *  1  is  calculated,  QJ  ♦  1  •  "  *  1  would  be  found 
from* 


QJ  ♦  1  ,  n  *  t  m  QJ  *  I  ,  »  +  AQJ  ♦«.»*»  7.12 

and  once  qj  ♦».*♦»  is  found,  the  real  solution  qJ *  1 • " *  1  would  be 
found  from  7.8  as:1 

qt  J  ‘  ‘  »  C w( Qu  ♦  i  ♦  Qu  -  t  )  /4  (l-w/2)Qu  1*  *  *  •  "♦  1 

7.  13 

When  solving  the  block  tridiagonal  form  of  7. 10  or  7. 11, 
the  solution  is  in  terms  of  the  vectors  AQ"  * 1 . 

The  solution  at  level  n,  which  is  used  to  calculate  the  right 
hand  side  of  7. 10,  is  in  terms  of  the  true  solution,  q.  The  only 
term  that  is  known  is  Afl" *  * ,  the  solution  of  the 
tridiagonal  matrix. 

AQ" ♦  *  must  be  changed  to 
Aq"*‘  before  it  can  be  added  to  q" .  In  other 
words,  AQ" *  1  is  smoothed  to  Aq" *  1 

which  is  then  added  to  q"  to  obtain  the  upgraded  solution 

q"*‘  *  q"  Aq"*1  7.14 
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This  is  vsry  similar  to  7. 12  but  with  Q  in  place  of  q.  If 
aquation  7. 13  at  laval  n*  1  is  subtractad  from  tha  same  aquation 
at  laval  n,  tha  following  is  obtained. 


A i  *  * .  »  ♦  *  ■  ♦  ^fli.  .  i  )•*♦*♦"*  1 /4 

♦  <  1-<V2>  AQl  J  *  ‘  •  "  *  ‘ 


7.  15 


4q"  *  1  is  now  in  terms  of  AQ*"  *  1 

which  is  tha  actual  result  of  tha  tridiagonal  solution  for  use  in 
Equation  7. 14. 

Equation  7. IS  is  used  from  nodes  L  ■  2  to  LMAX-1.  Care 
is  needed  at  nodes  2  and  LMAX-1.  At  node  LMAX-1 , 

AQU • t  is  £0  at  the  shock.  At 
node  2,  AQL - »  is  d3  at  the  wall 
and  both  must  be  accounted  for. 

At  node  2,  as  was  shown  in  Chapter  6: 


/Ap  t  \ 

1  AC  pT ) «  /T t  \ 

A(pu)  1 1 

0 

A<pw)tj  * 

G 

A<  pT) »/ 

A<  pT ) « 

w»  / 

l  Ap.  / 

At  node  2  then, 


Equation  7.15  becomes:  (For  components  4,  5) 


AqtJ  **•»**  ■  C<^flt*»/4  ♦  <  l-w/4)  AQl  ]J  ♦  1  •  "  ♦  1 


7.  IS 


For  component  1,  £Q».  - »  is  the  fourth  component  of 
£Q« ,  divided  by  wall  temperature.  For  components 
2  and  3,  7. 15  remains  the  same. 

At  node  LMAX-i,  as  was  also  presented  in  Chapter  6: 

■  Ql  ««*"**  “  Ql.lt  •*"  7.17 

This  is  known  at  the  shock  and  substituted  into  7.  15  during 
smoothing. 

The  calculated  4Q" *  ‘  in  the  tridiagonal 
algorithm  becomes  the  basis  of  convergence.  If  all 
«  *  approach  zero,  the  solution  at  the  current 
j+1  node  is  converging  properly.  The  criteria  used  here  was  that 
all  ABi*  *  *  must  be  less  than  or  equal  to  0. 000J 
before  moving  on  to  the  next  streamwise  node. 


8.0 


SXABmTX^NP.  EISHjVAfcUg  ANALYSIS 


Since  this  method  la  implicit,  it  is  expected  that  for  | 

reasonable  step  sizes:  fl 

Ax  or  Af 

the  finite  difference  algorithm  will  be  unconditionally  etabie. 

In  Reference  1,  the  step  size  variations  were  said  to  be  related 

to  changes  in  grids,  shock  propagation  accuracy  and  solution 

convergence  rates1 .  Ho  detailed  algorithm  was  given.  However, 

the  stability  considerations  to  be  discussed  here  are  not  related 

to  the  discrete  differencing  but  to  the  actual  physical  nature  of 

the  problem,  i. e. ,  has  the  current  formulation  eliminated  the 

elliptic  (subsonic)  region  of  the  domain?  ! 


Recall  Equation  5.7: 


CF,  J  *  1  •  -  ♦  At  *  ♦  1  •  •  Aq"*  1  3,  $ 

♦  t  F«  *  *  1  •  "  ♦  A«  *  *  1  •  •  Aq"*  1  3,  £ 
»€CSJ  *  1 «  "  ♦  MJ  *  1  •  "  Aq"  *  1  3,  £ 

♦  CH*  *  1 • *  ♦  A0J *  * • •  Aq"*  1  3 


8.  1 


Where  A» ,  Aa ,  M,  and  Ao  can  theoretically  change  in  both  the 
t  and  £  directions.  Now  assume,  as  in 

Reference  1,  that  these  coefficients  are  frozen.  This  is 
reasonable  if  changes  from  n  to  n+i  are  not  large.  The  above 
equation  is  now: 


i 
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-€m£q,  £  -  a0  Aq  •  G 


whvrt  at>  *i ,  m,  and  a0  art  frozvn. 


a.  2 


Am  m  further  simplification, 1  examine  the  viscid  and 
inviscid  limits  saparatsly.  Equation  8.2  is  too  difficult  to 
analyze  as  is.  Ths  inviscid  and  viscous  limits  art  simpler  and 
may  be  analyzed  separately.* 

8.1  INVISCID  LIMIT 

The  inviscid  limit  of  8. 2  can  be  written  as 


at  £q,  t  *  a«£q,  4  ♦  K<£q,  f,  £>  ■  0 

The  stability  analysis  of  the  system  of  equations  now  requires 
that  the  inverse  of  a*  be  formed  and  multiplied  throughout 
giving: 


Aq.  £  ♦  St'‘a*£q,  £  ♦  at  “  ‘  K  ■  0  8.3 

The  eigenvalues  of  the  a«.  ~ ‘  a*  matrix  now  determine  the  marching 
stability  of  equation  8. 3. *  Even  though  8.3  is  a  much  simplified 
version  of  5. 9,  it  is  a  form  suitable  to  mathematical  analysis. 

If  8.3  is  stable  then  the  full  equation  may  also  be  stable;  if 
8.3  is  unstable,  the  full  equation  certainly  will  be  unstable. 


As  •  further  simplification,  ihum  *  rectangular  grid 


for  this  analysis  eo  that  the  metrics  become 

fit  -  1  ft  ■  0 
4-0  4  ■  1 
J  »  1 


The  state  vector  is  unahang'ed. 

Introduce  the  following  differential  equation  of  state:4 

yp  -  ffT  +  0<P,  f  ♦  P,  4)  ■  0  a.  4 

This  becomes  a  generalised  fifth  equation  in  the  system  as  was 
discussed  in  Chapter  3.  Without  8. 4  the  Ft  and  F«  vectors  become 

pyi*  «■  p 

pvu  ♦  p 

<T/<r~l>"-V«  /2 ) pu 
0  / 
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F. 


r  \ 

4UV  ♦  P 
*  P 

(T/ ( y  1 )  +V*  /2>(0» 


The  matrioas  a>  and  a«  art  formed  from  the  partial  darivativaa  of 
tha  flux  vactora  with  respect  to  tha  atata  vactor  (Appandix  B). 

a 

An  eigenvalue  datarminatlon  raquiraa  tha  lnvaraa  of  at 
and  multiplication  by  that  invaraa.  Hovavar,  tha  laat  row  of  a« 
eonaiata  of  zaroa  ainca  tha  laat  componant  of  F»  ia  zaro  ao  tha 
invaraa  of  at  ia  undafinad.  To  ciroumvant  thia,  Bhutta  and  Lewis 
invantad  Equation  8. 4. 

Uaing  8. 4,  tha  fifth  oomponant  of  F»  and  Ft  bacomaa  6p 
(Equationa  3.16)  and  tha  laat  row  of  a«  and  a«  bacomaa 

0  0  0  0  6 


Now  the  invaraa  of  «t  can  ba  formed  and  tha  eigenvalues  of  ur* a* 
can  ba  found.  Although  Reference  1  indicates  that  8. 4  vaa  not 
used  in  tha  actual  computes'  coda,  a  subsequent  paper  by  Bhutta 
and  Lewis1  1  does  show  solutions  for  different  valuaa  of  6  and 
seams  to  indicate  aquation  8. 4  waa  included.  This  explains  tha 
inclusion  of  this  aquation  of  atata.  The  intent  vaa  to 
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investigate  different  values  of  6  and  its  possible  effect  on  the 
solution* 


The  eigenvalues  of  at " 1 a«  (based  on  0>O  but  for  the  limit  0 
approaching  0  >  are  t 

1,  w/u,  v/u,  w/u,  v/u  8.  S  | 

Reference  1  incorrectly  indicates  zero  in  place  of  unity.  The  | 

eigenvalues  of  the  4x4  matrices,  i. e.  without  the  equation  of  1 

state  as  the  fifth  equation  of  the  formal  system,  are:*  | 

w/u,  w/u,  uw  *  a<u“  ♦w*  -  a*>I/*/(u*  -  a*)  8.6 

The  essential  difference,  of  course,  is  that  in  8.  S  the  I 

eigenvalues  are  always  real  so  that  hyperbolia/marching  behavior  1 

is  indicated.  In  subsonic  regions,  the  second  set  of  eigenvalues 
(Eqn.  8.6)  become  imaginary  and  elliptic  behavior  is  indicated. 

Assuming  for  the  moment  the  analysis  of  Reference  1  is 
correct,  the  new  scheme  appears  to  be  stable  for  marching.  In 
contrast  to  previous  PNS  schemes  which  make  pressure  assumptions 
in  the  subsonic  sublayer  to  eliminate  the  instability,  the  new 
method  requires  none.  In  other  words,  even  though  the  same  j 

physical  problem  is  solved  by  both  methods,  mathematical  j 

reformulation  seems  to  grant  stability  to  the  new  scheme  without 


v* U*  N\ M VTUf  "Jin*  XV  JTsJ<Vjr4XV4VYV >f 


any  assumptions  in  ths  0*0. 0  limit,  for  which  ths  8. 8  aiganvaluas 
ars  to  bs  sxpsotsd. 

a.2  viscqus  mm 

In  ths  viscous  limit,  aquation  8.2  baoomas: 

a«  £q,  t  -  €mAq,  {  *  K  *  0  8.  7 

Using  ths  sama  6  rationals  as  in  ths  inviscid  limit,  ths  viscous 
sigsnvaluss  ars1  < 

0,  0,  0,  3€p/CPr/ou),  10€p/(3pu>  8.8 

In  prsvious  PNS  schsmss  ths  viscous  sigsnvaluss 
ars* 


0,  4p/(3pu),  <r j ,  <r 4 

Where 

<r,  .  *«  p/(2p(u*-a*  )  )  { <u+D/<uE>  )±C  <u>D/(uE)  )•  -4D/E3‘  '  •  > 

D»>u*-aa  E»Prp/k 

and:  k  ■  cosfficisnt  of  thermal  conductivity. 

Ths  discussion  concerning  marching  stability  is  rimilar  tc 
that  for  ths  inviscid  limit.  Ths  only  differsncs  is  that 
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positive  viscous  tigsnvaiues  require  positive  u-velocity 
components.  This  says  that  reversed  flow  is  unstable,  in 
agreement  with  accepted  practice. 

8.3  DISCUSSION  AND  ANALYSIS 

To  this  point  the  results  and  eigenvalues  may  seem 
surprising  given  that  mathematical  reformulation  seems  to  remove 
the  marching  instability  when  8*0.  The  analysis  was  first 
presented  in  Reference  1  and  the  results  of  equation  8. 5  check 
mathematically.  Whether  or  not  such  an  analysis  was  justified 
will  be  examined  here  using  a  small  perturbation  and  linear 
stability  analysis. 

Stability  for  the  mixed  differential/algebraic  system  of 
Equation  8. 3  is  unclear  since  3^0  implies  that  the  necessary 
matrix  operations  for  eigenvalues  are  undefined.  When  8  is 
included  to  form  a  totally  differential  system,  the  8  approaching 
0  limit  dees  not  yield  the  correct  eigenvalues  (Eqn  8.6).  There 
is  a  paradox  involving  the  parameter  8,  presumably  related  to  the 
singular  perturbation  form  of  the  modified  state  equation.  This 
reasoning  leads  to  the  following  analysis. 


6.  11 


The  eigenvalues  of  A  are: 

X. ,  .  1± < 2-a*  >*'* 

and  are  Imaginary  if  a  >  2. 

The  eigenvectors  are: 

e,  *  Cl,  ,-Xt  3T  ,  e*  *  Cl,  -x*  3T 

If  €  la  not  zero,  the  eigenvalues  of  8. 9  become 

Xi,  ...  =  0,  1±<2>*  '•  8.  12 

and  the  eigenvectors  are: 

ei  *  C 1,  0,  1 3T  ,  e»  *C  1,  -x*  ,  0  3T  ,  e*  a C 1 ,  —  Xi  >  0  ] T 

The  eigenvalues  are  real  and  have  no  "time  like"  constraint. 

\ 

Similarly,  the  eigenvalues  of  Equation  8. 5  were  unconditionally 
real  after  fthutta  and  Lewis  used  their  9  not  equal  to  zero 
assumption  to  allow  matrix  inversion  and  took  the  limit  as  0 
approached  Go  0. 

The  preceding  discussion  confirms  the  Bhutta  and  Lewis 
algebra  leading  to  8. 5.  The  essential  point  is  that  the 
eigenvalues  appear  to  be  independent  of  6  or  0  even  though  €  or  0 
cannot  apriori  be  set  equal  to  zero  to  find  the  eigenvalues. 

What  are  the  proper  eigenvalues  as  €  or  6  approach  zero  in  view 
of  the  differing  results  of  8.11  and  8.12  (or  8.5  and  8.6)?  Is 
it  not  necessarily  correct  on  physical  grounds  bo  assume  the 
eigenvalues  when  €  *  0  and  €  approaches  zero  must  be  identical? 


Consider  the  implications  of  assuming  a  regular 
perturbation  based  on  €.  First,  differentiate  the  second 
equation  of  8.9  v. r. t.  x  and  the  third  v.  r.  t  y»  Eliminate  w*  y  to 
get  s 

U«  -  Vy  *  0 

Vi*  -  Uy  *  ♦  2Vy  y  ♦  <  a*  Uy  Vy  > /€  *  0 


Then  with  vY  from  the  second  equation  and  uY  *  from  the  first 
equation  of  8.9: 

U«  -  Vy  *0 

€tv**  -  Vy  y  ♦  2vv  *  ]  ♦  a*  Uy  ♦  v*  -uV  ♦2Vy  ■  O  8.13 

Assume  a  regular  perturbation,  say: 

\ . 

u  *  u,  *  €ui  ♦  6*  u«  ... 

v  *  Vo  ♦  €v»  ♦  €*  v«  ♦  .  •  .  8.  14 

and  examine  the  eigenvalue  basis  for  different  orders  of  €. 


Substituting  8.  14  into  8. 13  and  collecting  terms 
according  to  like  powers  of  €  gives: 


0<<=e  ) 


Ud  x  -  Vo  V  *  0 


8.  15a 


Vo  -X  -  Uij  y  2v«j  V  +  4*U0y  a  0 
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0(€‘  )  Uu  -  Viy  *  0 


v*  %  ~  U|  y  .  ♦  2v,  y  ♦  »*  Ui  »  '  *  -  ( v#  »  t  -Vo  y  »  t2vfl  i  *  )  8.15b 


0(6*  )  Ui«  vtv  »  0 


v«  I  -  u*  r  r  ♦  a*u«v  *-(Vt  n“V,,,+2Vu»  )  8.  15c 


and  so  on. 


Thus,  In  general,  for  0<6"  >,  each  set  of  equations  can  be 


written  as : 


L(qn  )  J  R  <  Vi,  -  i  ) 

where  the  L  operator  is  Identical  for  the  equations  of  all 
orders. 


The  point  is  that  if  the  introduction  of  €  allows  a 
regular  perturbation  as  in  8. 14,  the  eigenvalues  should  be 
identical  for  all  orders  of  8,  including  the  zeroeth  order.  The 
Bhutta  and  Lewis  eigenvalues  are  not  identical  for  either  the 
model  (€  *  0, «1)  problem  8.9  or  the  actual  (0  *  O, «1 )  problem 
5.9. 


It  must  be  concluded  that  the  Bhutta  and  Lewis 
eigenvalues  are  inconsistent.  The  normal  subsonic  sublayer 


a.  3.2 


Bhutta  and  Lewis  erred  in  their  inviscid  stability 

i 

analysis  since  in  Equation  8. 3  they  did  not  consider  the  term 
K(£q,  ?,£),  which  is  also  a 

function  of  the  state  vector,  must  also  be  included  in  the 
stability  analysis. 


Rewrite  equation  8.3  as 
Aq„  ♦  Bq*  ■  Cq 

where  q  9  Lp,  u,  w,  T,  p]T  .  Using  a  linear  Fourier 
stability  technique,  set 


Then 


q  a  qe‘ * * *  *  *  *  * 


det<kA  ♦  IB  -iC>  *  0 


Here  the  A  and  B  matrices  are  the  partial  derivatives  of  the  X 


and  Z  flux  vectors  with  respect  to  the  state  vector  above 
(Appendix  B).  The  C  matrix  is  the  partial  derivative  of  the 
state  equation  vector,  H,  with  respect  to  q  above. 


Suggested  by  Dr.  Michael  B.  Giles 


On  evaluating  the  determinant  it  is  found  that 

either 


an 


k  «  <i/0>  (u*  -a*  )/ua 
which  is  unstable  if  M<1 

k  »  Cuwl  ±  <<uvl>«  -  <s»  -  u«  >»  1«  •  *  1/  (a#  -u«  ) 

1  *  Cuv  ♦  a<  u*  ♦  v*  -  a*  >  3/<a*  -a*  ) 

which  are  unstable  if  M<1 


The  conclusion  from  both  this  and  the  previous  section 
fl. 3. 1  is  that  the  Bhutta  and  Lewis  method  does  retain  the  weak 
elliptic  region  in  the  subsonic  sublayer.  Unless  appropriate 
measures  are  taken,  their  method  should  break  down  at  some 


downstream  location. 
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RESULTS 


Based  an  the  conclusions  of  Chapter  8,  the  present  scheme 
should  prove  to  be  unstable  at  some  point  downstream.  However* 
the  iterative  and  added  equation  aspects  of  the  scheme  make  it  so 
much  different  than  previous  PNS  formulations*  that 

the  instabiJ  ity  may  manifest  itself  in  different  ways. 

In  order  to  test  the  performance  of  the  present  code  over 
the  wedge  geometry*  cases  were  run  at  Mach  numbers  3  and  15  and 
Reynolds  numbers  ranging  from  approximately  4x10*  to  lxlQ7  at  the 
IDP.  Freestream  data  corresponding  to  80,000  feet  altitude  was 
used  when  dimensional  temperature  was  required  in  the  Sutherland 
law,  and  to  define  a  typical  reference  length  from  Reynolds 
number.  The  Mach  number,  Reynolds  number  and  Prandtl  number 
completely  define  the  air  flow,  irrespective  of  altitude  and  are 
the  only  relevant  inputs.  A  Prandtl  number  of  0. 72  was  used. 

The  results  of  the  numerical  tests  are  summarized  in 
Tables  9. 1,  9. 2,  and  9. 3.  These  tables  represent  a  progression 
in  grid  clustering.  The  results  in  Table  9. 1  were  carried  out  on 
an  evenly  spaced  50  point  grid  with  the  first  node  point  2.0%  of 
the  shock  layer  thickness  from  the  wall.  Table  9. 2  used  a 
clustered  50  point  grid  with  the  first  node  point  0. 6%  from  the 
wall.  Table  9. 3  used  a  clustered  100  point  grid  with  the  first 
node  point  0.3%  from  the  wall.  The  clustering  parameter,  s,  cf 


Chapter  4,  was  JL.OxlO-**1  ior  Table  9.1  end  3.0  for  Tables  9.2  and 
9.  3. 


The  progression  from  smooth  to  clustered  grids  was 
selected  so  as  to  obtain  general  performance  characteristics  of 
the  scheme.  Certain  cases  were  chosen  frum  Table  9. 1  and  rerun 
on  tighter  grids  to  obtain  the  results  of  Table  9. 2.  In  this 
manner,  similar  cases  were  comparable  with  respect  to  their 
performance  covering  genera'l  characteristics  over  the  entire 
shock  layer  to  very  specific  detail  within  the  boundary  layer  and 
in  particular  within  the  subsonic  sublayer. 
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TABLE  9. 1  SMOOTH- OHIO  RESULTS 


ALL 

VARIABLES 

RECORDED 

AT  INITIAL 

DATA  PLANE 

CASE 

MACH 

REl 

£»  L 

A» 

ii 

♦STEPS 

FINAL 

NO. 

Za  m 

X 

1 

3 

1.6x10s 

0.04 

0.03 

0.  013 

>150 

5.  5 

2 

3 

2.3x10* 

0.  10 

0.  03 

0.  021 

18 

1.  54 

3 

3 

4. 3x10* 

0.  21 

0.03 

0.  045 

10 

1.  30 

4 

2 

9.7x10* 

0.  05 

0.03 

0.017 

70 

3.  21 

5 

15 

3.8x10* 

0.  05 

0.03 

0.  0003 

>100 

4.  00 

6 

3 

1.6x10* 

0.  04 

0.001 

0.  013 

17 

1.  017 

7 

3 

1.6x10" 

0.  04 

0.  oni 

0.  013 

47,  9«0.01 

1.047 

8 

3 

1.6x10" 

0.04 

C .  001 

0.013 

60,  e»o. 10 

1.  069 

9 

3 

1.6x10" 

0.  04 

C.  06 

0.  013 

62 

5.92 

10 

15 

6.0x10* 

0.  04 

0.06 

0.  0003 

>100 

7.  00 

11 

3 

1.6x10" 

0.  04 

0. 10 

0.013 

50 

6.00 

12 

3 

1.  6x10* 

0.02 

0. 10 

0.013 

51 

6.  10 

13 

15 

1. lxlQT 

0.  03 

0.03 

0.  0002 

>160 

5.80 

14 

15 

8.9x10* 

0.  10 

0.03 

0.  0006 

15 

1.  4? 

15 

15 

1.  1x10* 

0.25 

0.03 

0.  00128 

3 

1.09 

16 

3 

1.6x10" 

0.  04 

0.03 

0.  006 

10 

2.  30 

17 

3 

2.3x10* 

0.  10 

0.  03 

0.019 

8 

1.  24 

18 

3 

4. 3x10* 

0.  21 

0,03 

0.05 

1 

1.  03 

19 

3 

4.8x10* 

0 ,  07 

0.  01 

0.  0125 

30 

1.  30 

20 

3 

1.  5x10" 

0.  04 

0.01 

0.006 

30 

1.  30 

21 

3 

2.5x10* 

0.  01 

0.01 

0.  004 

>100 

2.  00 

22 

15 

6.  0x10* 

0.  04 

0.  001 

0.  00017 

12 

1.  012 

23 

15 

6. OxiO* 

0.  04 

0. 0001  0. 00017 

>55 

1. 0055 

TABLE  9. 3  1QG  POINT  CLUSTERED  GRID  RESULT 
24  IS  6.0x10*  0.04  0.0001  0.00015  55  1.0055 
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Am  the  grids  become  increasingly  clustered  near  the  wall, 
the  apparent  height  of  the  subsonic  layer  changes  slightly  lor 
equivalent  boundary  layer  thicknesses.  This  is  due,  of  course, 
to  corresponding  node  location  changes  near  the  wall  as  the 
number  of  nodes  increases.  Cases  1,  5,  10,  13,  21,  and  23  all 
could  have  continued  far  a  larger  number  of  stops  than  indicated. 
That  does  not  mean  that  an  eventual  break  down  is  avoidable.  For 
example,  case  groups  <l,  2  and  3)  and  (13,  14  and  IS)  show  the 
resulting  behavior  for  increasing  boundary  layer  thickness  with  a 
constant  step  size.  Cases  1  and  13  were  concluded  at  a 
sufficiently  large  X  such  that  the  information  gained  showed  a 
trend. 

Figure  9. 1  shows  a  plot  of  the  number  of  marching  steps 
before  breaking  down  versus  the  height  of  the  initial  subsonic 
layer  for  Mach  *  3  and  a  stepsize  of  0. 03.  Figure  9. 2  is  a 
similar  plot  for  Mach  number  15  and  stepsize  of  0. 03.  The  same 
results  are  evident  in  cases  16  through  21.  A  dramatic  riae 
occurs  in  the  number  of  possible  steps  at  constant  stepsize  as 
the  subsonic  layer  decreases.  The  exponential  like  behavior  as 
the  subsonic  layer  shrinks  suggests  an  unlimited  number  of  steps 
would  be  achievable  for  an  inviacid  flew. 

Cases  9  and  10  provide  related  information.  All  other 
things  being  equal,  the  aoluticn  can  be  marched  farther 
downstream  at  high  Mach  numbers  because  the  subsonic  layer  is 
smaller.  This  is  apparent  in  cases  4  and  5  as  well. 
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It  la  interacting  to  compare  the  eclutiona  for  Mach  3  end 
Hech  13  in  ceeee  1  end  13.  Figures  9.  3  end  9. 4  show  the  profilee 
of  etete  verieblee  end  Mach  number  ecroee  the  ahock  layer  at  150 
atepa  for  Mach  3  end  13  respectively .  Notice  how  ah*»rp  the 
transition  ;.:one  between  the  boundary  layer  and  the  inviacid 
external  layer  has  become  when  compared  to  the  initial  data 
(Figs.  6.1  and  6.2). 

Figures  9. 5  and  9. 6  show  the  right  hand  aide  of  the 
governing  equations  (Cq.  5.6)  versus  vertical  node.  In  both 
cases  the  energy  equation  shows  the  largest  departure  in 
magnitude  from  zero,  followed  by  the  X- momentum  equation.  Recall 
that  the  right  hand  side  should  approach  zero  in  the  limit  of 
convergence  at  each  step.  A  nonvanishing  r^ght  hand  value  at  any 
node  indicates  the  extent  to  which  the  governing  equations  are 
not  being  satisfied  at  that  node.  Typically  it  was  found  that 
the  right  hand  side  was  of  0(1)  for  the  first  iteration  of  the 
first  streamwise  step  after  the  IDP  and  became  smaller  proceeding 
in  the  stramwise  direction  until  the  scheme  broke  down.  The 
largest  values  remain  within  the  boundary  layer*  and  are  caused 
primarily  by  an  inconsistent  set  of  initial  conditions. 

Figures  9. 7  and  9. 8  show  the  shock  surface  and  wall 
pressure  for  the  two  cases.  The  inviscid  shock  has  been  added 
for  reference.  The  large  jump  in  wall  pressure  could  be 
attributable  to  errors  in  the  initial  pressure.  Finally,  Figures 


9. 9  and  9. 10  are  enlargement*  of  the  boundary  layer  region*  for 
the  two  cases.  In  both  aaeee  the  reversed  density  gradient 
should  be  noted  and  in  Figure  9. 9  note  that  the 
gradient  is  lessening  and  eventually  reversed  flov  must  be 
anticioated. 

In  foot,  it  vai  found  that  independent  of  the  break  down 
location,  the  sheme  usually  broke  down  in  a  similar  way.  Figures 
9. 11  through  9. 17  trace  the  development  of  case  4  for  selected 
downstream  locations  (X  ■  1.12,  1.30,  1.90,  2.80,  3.10).  Notice 
particularly  the  behavior  of  pu  and  p. 

At  step  30,  Figure  9. 13,  the  density  has  a  sharp  gradient  and  the 
pu  variable  has  less  of  a  gradient  near  the  wall. 

Steps  SO  and  70  show  the  u-veloaity  approaching  separated  flow. 
Figures  9.1?  and  9.17  show  the  right  hand  side  values  <Equ.  3.8) 
for  viteps  60  and  70  after  they  have  converged.  Notice  that  the 
magnitudes  in  Figure  9. IS  show  that  the  governing  equations  are 
being  satisfied  while  those  in  Figure  9. 17  indicate  that  the 
solution  is  beginning  to  break  down.  In  the  same  manner  all 
cases  but  ona  which  become  unstable  do  so  ny  converging  to 
separated  flov.  Figure  9. 17a  details  the  boundary  layer  at  the 
onset  of  reversed  flow  for  Case  20.  The  large  density  spike  Is 
characteristic  of  the  solution  oppraachlng  separation.  The 
governing  equations  remain  satisfied  until  the  break  down  is 


imminent 


Figure  9. 1C  shows  ths  distribution  of  pressure  across  the 
shock  isysr  for  case  9,  which  also  canvcrgu  to  separated  flow. 
Ths  profiles  art  plotted  at  stsps  2,  4,  8,  16,  32  and  64  down  ths 
wsdgs  (X-1.12,  1.24,  1.48,  1.96,  2.92,  4.64).  Figurs  9.18a  is  a 

blowup  of  9.18.  Ths  profils  at  stap  64  (X  *»  4.84)  is  approaching 
ths  prsasurs  profils  of  Figurs  9. 3  which  had  gons  13G  stsps  with 
a  slightly  higher  Rsynolds  numbsr.  This  stems  to  show  that  for 
nsarly  idsntical  initial  conditions,  ths  scheme  prsdicta 
canaiatsnt  solatio.  j  independent  of  ths  two  different  streamwiae 

A 

stspaizss  used.  (A^*0.Q3  in  Figure  9.3  and 
Ax*0.06  in  Figurs  9.18). 

i 

Anothsr  interesting  phenomenon  ic  observed  with  an 
unstable  Ax.  Consider  cass  6.  This  is  the  only 
case  that  did  not  break  down  due  to  flow  separation.  Figures 
5.  19  and  9. 19a  show  that  pu  is  tot  tending  toward  a 
decreasing  gradient.  It  in  also  interesting  to  ate  that  the 
errors  are  at  or  below  the  sonic  line,  which  is  shown  for 
reference.  A  possible  reason  for  this  is  suggested  in  Chapter 

10.  Figure  9.20  shows  the  growing  magnitudes  of  the  right  hand 

£ 

aide  of  the  governing  equations  (Cq.  5. b)  for  case  6. 

Finally,  consider  case  24  of  Table  9. 3.  Recall  that  this 
result  was  obtained  on  the  most  clustered,  100  point,  grid. 
Previous  PNS  schemes*  •  1  * *  ‘ 8  have  documented  the  phenomenon  of 
departure  solutions.  Departure  behavior  is  found  in  PNS  codes 
that  do  not  include  a  pressure  assumption  and  is  characterized  by 


an  exponential  drop  in  wall  pressure  associated  vith  the  subsonic 
ellipticity.  Figure  9.21  shows  the  wall  pressure  distribution 
lor  case  24  and  Figure  9.22  shows  a  u-velocity  component  profile 
for  two  downstream  X  locations.  Both  figures  are  consistent  with 
previously  demonstrated  departure  behavior* • *  * • *  8  including  the 
rapid  development  of  the  pressure  drop.  Figure  9. 23  shows  the 
right  hand  side  distribution  from  equation  5. 8  and  indicates  that 
the  largest  errors  appear  in  the  energy  equation.  This  again  is 
consistent  with  departure  behavior*  •  ‘ *•  ‘ a .  Figure  9.24  displays 
the  lower  portion  of  the  shock  layer  at  departure.  Note  the 
large  negative  gradient  in  /ow.  Compare  this  with 
Figures  9. 23  and  9. 28,  which  are  from  case  23.  The  only 
difference  betveen  the  two  cases  is  the  grid  clustering.  They 
are  displayed  at  the  same  downstream  location,  X=1.0055.  Despite 
this,  case  24  indicates  departure  while  case  23  does  not.  Note 
particularly  that  the  negative  pv  gradient  exists  in 
figure  9. 25.  But  the  pressure  and  energy  behave  quite 
differently  than  those  in  Figure  9. 24.  This  is  consistent  with 
and  ap - ?ars  to  validate  the  claim  of  Reference  15  which  states 
chi.it  mesh  spacing  and  stepsize  must  be  of  0<Re“3'*)  or  less  in 
order  to  pick  up  departure  behavior. 

The  indications  from  cases  x-5  and  9-23  is  that  even  when 
classic  departure  behavior  is  not  evident,  the  instability 
manifests  itself  by  converging  to  e  jarated  f lo".  Figure  9.27 
snow*  the  right  hand  side  distribution  from  equation  5.8  ox  case 
23.  Of  importance  is  that  even  though  departure  behavior  is 

8i 


suppressed  errors  do  grow  in  the  energy  equation,  which  is 
consistent  with  case  24. 

Figure  9. 23  provides  a  last  comparison.  Pressure 
distribution  profiles  in  the  lower  part  of  the  shock  layer  are 
shown  at  different  X  locations  during  departure.  Notice  that  as 
the  solution  marches  downstream  a  pressure  disturbance  moves 
outward  away  from  the  surface  in  the  shock  layer.  This  may  be 
compared  with  the  behavior  in  Figure  9. 17,  which  was  run  at  Hach 
3  and  also  exhibits  an  outward  moving  pressure  wave.  Figure  9. 17 
represents  a  somewhat  more  stable  solution  since  as  the  solution 
marches  downstream,  the  pressure  tends  towards  constancy  and  is 
consistent  with  known  wedge  flow.  Pressure  in  Figure  9. 23,  on 
the  other  hand,  moves  away  from  constancy  but  the  underlying 
instability  in  both  cases  possibly  induces  similar  behavior. 

The  cases  shown  have  been  selected  because  they  provide  a 
representative  characterization  of  the  behavior  of  the  current 
scheme.  They  appear  to  imply  that  the  scheme  is  unstable,  that 
it  is  unstable  due  to  departure  behavior,  and  even  when  classic 
departure  behavior  is  not  evident  due  to  grid  sizing  or  step 
size,  the  underlying  instability  causes  convergence  to  separated 
flow. 


The  results  are  not  exhaustive  but  they  do  serve  to 
outline  the  performance  of  the  current  algorithm.  Because  of  the 
iterative  nature  of  this  particular  algorithm  (Chap.  3)  the 


classical  PNS  pressure  modeling  could  not  be  included.  Several 
of  the  previous  sublayer  assumptions  were  mentioned  in  Chapter  1 
and  it  would  have  been  informative  to  use  similar  sublayer 
assumptions  in  the  current  scheme,  hoping  to  achieve  a  stable 
algorithm. 
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Figure  9.3  State  Variable  Profiles  after  150  Steps  at  Mach  3 
Ax  ■  0.03 
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Figure  9.5  Values  of  Governing  Equations  after  150  Steps  at 
Mach  3,  Ax  -  0.03 
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Figure  9.7  Wall  Pressure  and  Shock  Surface  after  150  Steps  at 
Mach  3,  Ax  =  0.03 
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Figure  9.10  Boundary  Layer  after  150  Steps  at  Mach  15,  Ax  =  0.03 
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Figure  9.12  Case  4  Profiles,  X  *  1.30 
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Figure  9.15  Values  ox  Governing  Equations  of  Case  4,  X 
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Figure  9.17a  Onset  of  Separated  Flow,  Case  20 
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Figure  9.18  Streamwise  Pressure  Distribution  of  Case  9 
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Figure  9.18a  Blowup  of  Case  9  Streamwise  Pressure 
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Figure  9.20  Values  of  Governing  Equations  for  Case  6 
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Figure  9.22  Departure  Velocity  Profiles,  Case  24 
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Figure  9.23  Values  of  Governing  Equations  at  Departure,  Case  24 
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Figure  9.25  Boundary  Layer  Profiles,  Case  23 
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Figure  9.28  Streamwi9e  Pressure  During  Departure,  Case  24 
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10.0  CQHCLUSIQKS 


The  Dhutta-Lewis  Algorithm  (Ref.  1)  b«hav*i  in  an 
expeated  way  for  PNS  jchemew  that  do  not  make  apaoifie  aublayar 
aaauaptiona.  Of  particular  ralavanoa  to  thia  work  in  reaahing 
thia  conclusion  ara  Rafaranaaa  13,  K.  and  IS.  Tha  lattar 
apaoifioally  addraaa  tha  role  of  PNS  aodaling  and  ralatad 
quaationa  of  marching,  dapartura  bahavior  and  atability  for  PNS 
aquationa.  Thay  particularly  diacuaa  in  aoaa  dapth  tha  usaa  of 
tha  atraaawiaa  atap  aiza  to  ovarooaa  waak  allipticity.  Tha 
praaant  analyaia  uaaa  atap  aiza  aa  a  way  to  eharaatariza  tha 
bahavior  of  tha  Bhutta-Lawia1  approaoh. 

Tho  raaulta  of  Chaptar  9  indicata  that  tha  algoritha 
convargaa  to  and  breaks  down  by  approaching  a  separated  flow 
condition.  Thia  is  in  agraaaant  with  tha  results  of  tha 
stability  analyses  of  Chaptar  8.  Exponential  dapartura  bahavior 
may  not  ha  apparent  unless  tha  grid  and  atepaize  ara  sufficiently 
•mall,  but  tha  dapartura  bahavior  in  manifested  by  convergence  to 
separated  flow. 

Recall  from  Chaptar  S  that  tha  tridiagonal  algorithm 
procedure  ia  associated  with  pauedo- temporal  steps  (iterations) 
rather  than  apace.  Therefore  tha  atability  analyses  of 
Rafarancas  13,  14  and  IS  which  involve  implicit  space  steps  ara 
not  applicable.  However,  it  ia  reasonable  to  assume  that  tha 
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ourrent  iterative  approach  aaounta  to  taking  the  aaae  atep  size 
repeatedly  and  atep  eize  doee  then  eharaoterize  the  behavior. 

A  ainiaua  atreaava.ee  etep  eize  ia  related  to  the  height 
of  the  aubaonic  aublayer . * •  •  1  * • *  *  The  relation  developed  in 
Reference  13  1st 

iX/i|  >  i/R 

where  &•  ia  the  thiokneaa  of  the  aubaonio  layer.  The  larger  the 
elliptio  zone,  the  sore  restrictive  on  aoouraoy  the  above 
oondition  becoaea. *  *  Thia  explaina  the  reaulta  of  Caae  6.  Also, 
it  ia  iaplioit  that  atable  aolutiona  are  aore  easily  achieved 
with  high  speed  flows  in  view  of  their  relatively  smaller 
aubaonic  layers. 

This  behavior  is  consistent  with  the  findings  described 
in  Chapter  9.  The  algorithm  will  aaroh  downstream  when  certain 
conditions  are  met.  Only  if  the  subsonic  layer  is  sufficiently 
small,  e.  g.  as  in  the  4X  or  5X  boundary  layer  cases  of  Tables  9. 1 
and  9.2,  is  stability  evident.  Larger  boundary  layers  can  be 
aarchod  with  higher  Mach  numbers  because  the  subsonic  layer 
remains  small.  This  agrees  with  References  13  and  14.  In  fact, 
when  this  condition  is  met,  a  rather  large  step  size  of  0. 06  was 
used  successfully.  Smaller  step  sizes  of  order  10*  *  or  10* 4  also 
work  and  step  sizes  of  this  order  have  been  used  by 
others. * • » 4 • * »  It  must  be  concluded  that  the  Bhutta-Lewis 
algorithm  (Ref.  1)  behaves  in  the  same  manner  as  was  to  be 
expected  for  PNS  schemes  without  sublayer  assumptions. 
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It  should  bo  noted  that  the  hypersonic  similarity 


parameter 


H./<Re>‘ '  ■ 


was  the  basis  here  for  selecting  specific  Reynolds  numbers 
(Chapter  9).  At  Hash  numbers  3  or  IS,  a  match  of  the  hypersonic 
similarity  parameter  to  that  of  Bhutta  and  Levis'  original  paper 
(Ref.  1)  implies  a  boundary  layer  thickness  exceeding  20X  of  the 
local  shock  layer  thickness.  Xecording  to  the  aformentioned 
stability  analysis, 1  *  this  should  severely  reduce  the  allowable 
step  size.  Indeed,  the  current  research  failed  to  march  low 
Reynolds  number,  viscous  flows  with  large  boundary  layers  any 
appreciable  distances  downstream.  Bhutta  and  Lewis,  however, 
developed  the  current  algorithm  to  solve  low  Reynolds  number 
flows  and  were  apparently  able  to  march  largr  distances 
downstream. 


One  portion  of  the  Bhutta  and  Lewis  algorithm  which  the 
current  work  did  not  duplicate  was  variable  step  size.  In 
Reference  1,  results  indicate  that  the  step  size  begins  with 
Q(10~a)  and  ends  with  0(1).  Other  documented  parameters  of  the 
Bhutta  and  Lewis  algorithm  have  been  duplicated  and  tested  and 
have  been  shown  to  not  ensure  marching  stability.  Since  the 
current  work  agrees  with  past  references  on  the  importance  of 
atreamwise  step  size,  it  may  be  that  variable  step  size  enabled 
the  instabilities  to  be  supreased.  This  is  consistent  with  the 
fact  that  successful  marching  was  indicated  with  high  Reynolds 
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number  flows  whsrs  ths  boundary  lsysr  and  subsonic  laysr 
thieknsssss  do  not  grow  appreciably  with  downstrsas  distanos.  As 
ths  boundary  laysr  and  subsonic  laysr  grow  in  a  viscous  flow,  ths 
stsp  size  apparently  must  increase  to  maintain  stability,  as 
implied  in  Reference  1. 

Nevertheless,  the  scheme  has  been  shown  to  be  unstable. 
Therefore,  in  order  tc  march  successfully,  some  parameter  or 
parameters  in  this  scheme  must*  be  adjusted  to  allow  for  marching 
stability.  Unfortunately,  the  parameters  which  allow  marching 
stability  were  not  discussed  in  Reference  1  and  no  reference  was 
made  to  the  characteristic  instability  of  the  scheme. 

The  role  of  the  differential  modification  of  the  equation 
of  state  is  also  interesting.  Recall  that  the  modified  equation 
of  statei 

HP  ♦  8<P, f  ♦  P, €  )  »pT 

was  used  to  replace  the  true  equation  of  state.  For  most 
calculations  8  was  set  identically  equal  to  zero;  in  certain 
cases,  8  was  allowed  to  vary.  In  Reference  1,  Bhutta  and  Lewis 
indicate  that  8  was  not  used  for  actual  numerical  cases,  but  in 
Reference  11  they  present  solutions  with  8»10"**  and  10~*  and 
conclude  that  such  magnitudes  have  negligible  effect.  The 
current  study  also  showed  no  influence  to  be  present  with  8 
chosen  to  be  a  small  number. 
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However,  for  larger  theta,  (oaeee  7  and  8)  say  8*0.01  or 
0. 1,  the  perforaanoe  vaa  actually  improved  in  the  aenae  that 
marching  oontinued  farther  to  larger  Xu**.  It  ahould  be  noted 
that  improvement  waa  noted  only  for  high  Reynolda  number  and 
email  boundary  layer  cases |  namely,  the  thin  aubaonio  aublayer 
case.  Appreciably  viacoua  flowa  were  unchanged. 

10.1  AREAS  FOR  FUTURE  RESEARCH 

Three  auggeationa  for  poaaiblo  further  reaearch  into  the 
algorithm  inalude  the  uae  of  a  aublayer  assumption,  a  nonuniform 
atepaize  and  more  aanaiatent  initial  conditiona.  The  acheme 
ahould  be  explored  with  appropriate  preaaure  aaaumptiona  inaerted 
during  the  derivation  of  the  block  tridiagonal  form.  Thia  will 
allow  inveatigation  of  the  performance  of  thia  acheme  with  the 
claaaical  aublayer  aaaumption. 

Aa  mentioned  earlier,  the  algorithm  of  Bhutta  and  Lewie 
uses  a  variable  atep  aize  baaed  on  aeveral  factora.  In  the  very 
viacoua  caaea  where  the  atable  atep  aize  regime  ia  email, 1 9 
perhape  a  correctly  calculated  atep  aize  would  provide  a  more 
aatiafactory  marching  performance.  Aa  the  boundary  layer  and 
aubaonic  layers  grow  in  the  atreamwiae  direction,  the  atep  aize 
can  grow  accordingly.  Thia  ia  evidenced  in  Reference  1,  but 
there  was  no  clearly  stated  methodology  to  indicate  how  the  step 


isf .ze  vaa  calculated  or  how  it  varied 


Figure  10. 1  shows  the  constraint  schematically.  As  the 


subsonic  layer  Increases  from  level  1  to  2,  the  minimum  step  size 
also  Increases.  Assuming  a  maximum  step  size  for  a  given  level 
of  accuracy,  there  Is  a  finite  step  size  window  which  decreases 
as  the  subsonic  height  Increases  until  there  Is  no  stable  step 
size  available. 

Obtaining  a  better  Initial  condition  could  expand  on  the 
results  of  this  study.  Bhutta  and  Lewis  used  a  starting  code 
whereas  this  work  uses  only  an  estimate  for  initial  conditions. 
Indications  found  during  computer  runs  for  Chapter  9  showed  that 
initial  conditions  which  better  satisfy  the  governing  equations 
may  allow  the  solution  tc  march  farther,  all  other  factors  being 
equal.  The  computer  runs  also  showed  that  the  initial  profiles 
for  low  and  high  Mach  numbers  are  not  necessarily  similar. 

The  initial  conditions  used  in  this  work  are  a  compromise 
between  Mach  numbers  3  and  15.  Different  initial  conditions  at 
each  Mach  number  would  have  satisfied  the  governing  equations 
more  completely  but  this  was  not  deemed  to  be  crucial.  As  an 
illustration,  return  to  Table  9.  1,  case  10.  The  initial 
conditions  (Fig.  6.2)  provided  for  a  slightly  nonlinear  profile 
slope  from  the  wall  to  the  top  of  the  boundary  layer.  When  this 
was  changed  to  a  linear  slope,  Case  10  was  run  for  200  steps  at 
Axa0. 06  to  a  final  X  position  of  13. 0.  The 

results  are  shown  in  Figures  10. 2  through  10. 5.  Judging  from  the 
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■mall  right  hand  aide  magnitudes  o f  Figure  10. 4,  the  solution  was 
still  satisfying  the  governing  equations  and  would  have  marched 
farther.  A  more  consistent  set  of  initial  conditions  would  not 
change  the  stability  characteristics  of  the  scheme  but  they  would 
remove  an  extra  source  of  error,  which  would  aid  in  marching  the 
solution  downstream. 


Figure  10,2  State  Variable  Profiles,  X 
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APPENDIX  A 


EQUATION  TRANSFORMATION  DERIVATION 
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Wi 


Host  of  the  basic  transformation  from  a  raotangular  to  a 
gsnaral  coordinate  system  vas  performed  in  Chapter  3.  The 
transformed  equations  are  3. 13a  and  3. 14.  The  final  step  is  to 
transform  the  individual  components  of  the  vectors  in  the 
governing  equations. 

Comparing  3. 4,  3.13a,  and  3.14  ve  see  that: 


F»  -  fid  ♦  ft  Gi 
Fa  -  faEi  ♦  6  61 


A.  1 


Substituting  for  Ei  and  Qi  above,  the  following  forma  result  for 
Ft  and  Fa  : 

/p<f«u»f«w>  \ 


F,  «1/J 


pu<f.  U*fa  vWiP 
pw<f,  u*|*  w>«-f«P 
<T/  <  7-1 )  /2  >p<  f«  u*fi  w ) 

i  0 


/ 


A.  2 


I 


F.  «  1/J 


p<& i  u-^it  w) 

pw<4t  u+t*  v>  ♦  {«  P 
( T/ <  7-1  >  ♦  V* /2 )  p<  fi,  u-^fii  w  > 
0 


\ 


\ 


A.  3 
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If  the  Modified  differential  state  equation  is  used,  the 
final  component  of  Ft  and  F«  will  be  6p.  The  1/J  factor  is  added 
to  preserve  the  strong  conservation  lav  form. *  It  is  differenced 
with  the  other  components  of  the  vector  which  preserves  the 
conservation  property.  Note  that  in  Reference  1  the 
contravariant  velocities  t 

<fk  u+f(  v>  and  (&u*&v) 
were  printed  incorrectly. 


The  last  vector  to  be  transformed  is  the  viscous  vector 
since  the  vector  H  is  algebraic  and  is  not  transformed.  The 
viscous  veator  is  slightly  more  complicated  than  Ft  and  F* 
because  it  contains  X  and  2  derivatives  which  must  also  be 
transformed.  From  Equations  3.4,  3.13a,  and  3.14  again* 

S  ■  (m  Ev  ♦  fii  Gv 


where  Gv  was  seen  in  equation  3. 4  and 


/ 


ti  i 

Ev  »  I  T*  1 

ut» «  ♦  VT«  1  -  q» 

0 

The  derivation  will  be  carried 'out  component  by  component. 
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Tha  first  eomponant  oi  S  is  zero.  Tha  aacond  eomponant 


1st 

Ci  *Ti  i  ♦  Cl  *T*  i 
vhieh  la  axpandad  sat* 

pCC*<4u,  i  -  2v,  t>/3  ♦  &(*,«  ♦  u,  t>3 
Thaaa  X, 2  darivativaa  muat  alao  ba  tranaf ormtd  to  the 
f  ,  C  ayataa. 

Using  3.12  tha  abova  axpanaion  ia  tranaf ormad  tot 

p(Ci  C4(fiU,  f  ♦  Ci  u,  0/3  -  A.  4 

2(fi«f!  ♦  Ci  w,  0/33  ♦ 

CiCfiW,  t  ♦  Ci  v,  C  ♦#*  u,  f  ♦&  u,  C3  > 

Following  tha  PNS  assumptions,  all  straamwiaa  viscous 
darivativaa  ara  droppad  from  A. 4  laavingt 

(p/J)CM0u,  C  ♦  (Hi  iu,  C  ♦  Mt»w,  C)/33  A.  5 

whara 

M.  -  &•  ♦  C*' 

M. ,  -  &• 

Hi  i  a  C& 

Tha  third  componant : 

Cl*Ti«  ♦  Cl  *Ti  i 

ia  dona  in  tha  aama  mannar,  rasulting  in: 

<p/J>CMoV,C  ♦  (Hi  i  u,  C  ♦  Hi  i  w,  £)  /33  A.  6 
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where 


where 


Mi  *  *  St  * 

The  fourth  component  let 

£■  <u*tii  ♦  WTt*  -  q*)^Ct(u'Tti  ♦  WTi »  -  qt  > 


q»  ■  -k(T,  i  )  q«  ■  -k(T,  »> 


The  final  form  of  the  entire  v‘iacoue  vector  ie  given  byi 


CM«u,  {  ♦  (Mi.u,  £  ♦  Mu  w,  0/3] 

CM*  w,  {  ♦  (M,  iU,  {  ♦  Hi  i  w,  {)  /32 

CM*  (T,  {/(Pr(yl>  >>uu,  {♦ww,  {) 

♦  (Mu  uur  {♦Mi  i  ww,  {♦Mi  «  (  wu,  {♦uw,  {>  )  /3 3 


In  Reference  1,  the  first  two  of  the  1/3  coefficients  of  S  were 
incorrectly  printed  as  1/2. 
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APPENDIX  B 


JACOBIAN  MATRICES 


After  the  coordinate  transformation,  F, ,  F« ,  S,  and  H  are  given 
by  3. 15  , 3. 16  and  Appendix  A.  The  state  vector  after  the  same 
transformation  is 

q  ■  1/J  Ep,  pu,  pw,  pT,  p]T  B.  1 

The  Jacobian  matrices  are  formed  by  the  following  partial 
derivatives: 


A«  *  H,  q  A,  ■  F»,q  A.  ■  F« ,  q  M  ■  S,  q  B.  2 

For  example,  the  first  row,  second  column  of  Matrix  A.  is  the 
partial  derivative  of  the  first  component  of  F»  with  respect  to 
the  second  component  of  q  (Eq.  B. 1>.  The  fourth  row,  fourth 
column  of  matrix  Aa  is  the  partial  derivative  of  the  fourth 
component  of  F*  with  respect  to  the  fourth  component  of  q  (Eq. 
B. 1).  Each  component  of  q  is  treated  as  one  term  when  taking 
derivatives.  They  are  not  broken  down  into  primitive  variables. 

Following  are  the  four  Jacobian  matrices  using  the  modified 


equation  of  state. 


A* 


J 


0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  -1  7 


At 


0 

f* 

ft 

0 

-uUt 

2fe  »fe  w 

ft  u 

0 

-wUt 

it  w 

ft  u+2ft  w 

0 

T, 

T«  ♦$*  uw 

T»  +  ft  uw 

Ut  /  ( 7~  1 ) 

0 

0 

0 

0 

where 

Ut  ■  feu  +  ft* 

T,  »  -U,  (T/  (7-1)  ♦  V«  ) 

T*  «  ft  <T/<7-1>*(3u*  ♦w*  )/2) 
T,  =  it  <T/<7-l>+<u*  +3w«  )/2> 


At 


r  ° 

ft 

ft 

0 

i 

c 

c 

» 

2ft  +ft  w 

ft  u 

0 

-wUt 

ft  w 

ft  u+2ft w 

0 

T  4 

T«  +ft  uw 

Tt  +ft  uw 

Ut  /(7-1) 

0 

0 

0 

0 

where 

Ut  ■  ft  u  ♦  it  w 

Tt  -  -U.  <  T /  <  7~  1 )  ♦  V«  ) 
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0 

ft 

ft 

O 

e 


o 

St 

St 

0 

e 


i 


i, 

w*M**Mwuwutniwwicuwwfiw\iw«ify»i  vww»rvi  junui  jwfjonrvxifi'  V^^i,7J(?J(^>XXAX.VUWJL"jA^r 


T.  ■ 

T.  «  &  <T/<7-l>*<u*«-3w*  )/2) 


M  »-p/J 


0 

BD|  *aDt 

<rD«  +aDi 

(ND«  ♦flDn 
♦  crD*  «-2aDT  > 
0 


0 

-flD* 

-aD* 

-flDj  -aD*" 
0 


0 

-aD* 

-aD* 

-aD*  -aD» 
0 


0 

0 

0 

-HD* 

0 


where 

<3  *  Mo+M*i/3 
a  *  Mx  *  /3 
a  «  Ha ♦Mx  x /3 
N  -  Mo  /  <  Pr  <  7- 1 )  ) 

The  following  are  all  partial  derivatives  contained  in  M: 
D,  -  (uj/is)t$ 

D.  -  (wJ /p),S 
D,  =  <J /P),S 


D*  -  (TJ /p),& 


D»  »  <u*J/p>,£ 
D.  »  ( w*  J/p) ,  fi 
D7  ■  ( uwJ/p) ,  £ 


0 

O 

0 

0 

O 
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APPENDIX  C 


MISCELLANEOUS  DERIVATIONS 
MASS  FLOW,  U, ,  W. ,  SHOCK  PREDICTION 


Equations  C. 1  and  C.  2  art  two  aquations  in  tha  two 


unknowns,  U«  and  Wa .  Solving  simultaneously  givasi 

AA  ■  1  ♦  tan*  <fl  -  n/2  -  6.  ♦  cos-  ‘  <  tanfl/  <  ft  •  ♦  tan*B>*'*> 


U.  »  (1/tan*  fl  ♦  l/ft*)/AA3‘'*  C.  3 


W.  -  U.  (AA)  C.  4 


MASS  FLOW 

Tha  geometric  relationships  far  mass  flow  are  shown  in 
Figure  C.  2.  The  freestream  normal  distance,  Z<a  ’ ,  must  be  found  in 
terms  of  the  wedge  Z  coordinate. 


Figure  C. 2  Mass  Flow  Geometry 
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For  mass  flow  to  balance.  It  must  satisfy: 


( puZ )  CD  » 


pi  ■  m  a  e  k 
J •  pudz 


After  nondlmeneionalizlng  and  finding  the  relationehip 


for  Zm,  Equation  C. S  becomea : 


M®(Z* m coe<  Gt  >  ♦  Xein(8i>>  ■ 


pudz  C.  6 


To  discretize  the  integral  in  C.  6,  the  vertical  nodes  <L) 


of  Chapter  S  are  used.  The  integral  is  then: 


El  <  pu  )  Az 


where  the  summation  is  carried  out  over  the  vertical  nodes. 


The  value  of  the  mass  flow  across  the  shock  (Equ.  C.7) 


normalized  with  respect  to  the  freestream  mass  flow  (Equ.  C. 6) 


should  be  near  1.0.  If  not,  the  shock  position  is  moved  up  or 


down  to  equalize  the  two  mass  flows.  Reference  1  suggests  that  a 


difference  of  0. IX  is  sufficient. 


SHOCK  PREDICTION 


From  the  reference  by  Chausee  et  al.T,  the  shock  location 


in  two  dimensions  can  be  predicted  from  location  j  to  j  +  1. 


Differences  in  coordinate  systems  have  been  accounted  for.  The 


m  she  location  is  given  by; 


Z*  *  ‘  -  ZJ  *  Z ,  f  Af 


where 


Z,  f  -  <T?y£i  -  17*  &  >/J 


which  in  tvo  dimensions  reduces  tot 


Mi  *  /  ( M®*  -  M,  „•>*'«  *  Mi  n  /Mi  «  *  tanfl  C.  9 


(Fig  C.  1).  Then  Equation  C.  8  bMQMi 


2*  ♦  *  ■  ZJ  ♦  tanflnx 


This  is  shown  geometrically  in  Figure  C.  3 


Figure  C. 3  Shock  Prediction 
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m  ■  m  « 


•«  w-w  w-m  kts  wm 


w  wivaum  wti  \fm  u  « m  ura  u*  L'v  u%  '.«im  vw.»y  lvi\'  lvw\  ^*\  ram  mx  vx  *A*.v1 


APPENDIX  D 


COMPUTER  CODE 


Figure  D.l  Flowchart 
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Th*  cod*  la  tltlad  PARABNS. FOR  and  uaaa  th#  £ 11a 


PARABNS. INC.  Both  ara  written  In  FORTRAN  and  hava  baan  run  on  a 
Digital  VAX  11/760  and  a  HloroVAX. 


SUBROUTINES 


INIT  -  Takaa  Input  of  tha  boundary  layar  thleknaaa  and 
aakaa  aatlaataa  for  tha  Initial  profllaa.  Uaaa  aaaaflov  to  raflna 
Initial  aatlaata  for  ahook  angle.  Inviacld  ahook  angla  and 
Ranklna-Hugoniot  ralatlona  provlda  Initial  estimates  at  tha 
ahock. 


ORIDGEN  -  If  tha  time  laval  la  at  tha  flrat  Iteration, 
tha  routine  predict a  tha  j+l,n«l  at at a  vector  baaed  on  paat  at ate 
vactora.  Alao  predict*  nav  ahock  location.  If  tha  time  laval  la 
greater  than  one,  mass  flow  la  uaad  to  refine  tha  ahock  location. 
Interior  value*  of  danaity  are  uaad  to  predict  tha  ahock  danaity 
and  therefore  tha  ahock  angle.  Other  ahock  value*  are  found  from 
tha  Rankine-Hugoniot  relation*.  Metric*  and  X, Z  poaitiona  are 
computed. 


ASSEMBLE  -  Aaaambla*  the  individual  Jacobian  matrices 
into  block  tridiagonal  form.  It  alao  perform*  tha  tridiagonal 
boundary  condition*.  Two  subroutine*  are  called) 

A)  FINDG  -  Calculate*  th*  right  hand  side 
vector  (Equ.  3.8). 

B)  BUILDMATRIX  -  Fills  the  individual  Jacobian 
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matrices!  A0 ,  A,,  A*  and  N 


BSSOLV  -5xS  Block  tridiagonal  aolvar 

FINDSTATE  -  Smooth**  th*  £0" *  1  tarm*  and 
updataa  th*  *tat*  v*otor  with  th*  newly  smoothed  solution. 

OUTPUT  -  Creates  5  output  files t 

A)  State  Variable  Profiles 

B)  Convergence  History 

C)  Governing  Equation  Errors 

D)  Shock  Surfaoe  and  Wall  Pressure 

E)  Streamvise  Profiles  of  Selected  State 
Variables 

A  copy  of  the  aomputer  cod*  is  on  file  with  Professor  Jud^on 
R.  Baron.  Massachusetts  Institue  of  Technology.  Room 
33-217,  (617)  233-4329. 
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